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sunspots  is  studied,  in  order  to  assess  the  possibility  of  cooling  by  Alfv£n 
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shows  strong  upward  reflection  In  the  convection  zone  as  well  as  strong 
downward  reflection  in  the  photosphere  and  low  chromosphere.  These  results 
tend  to  rule  out  significant  sunspot  cooling  by  waves. 

A  study  of  a  simple  thermal  model  of  a  sunspot,  based  on  the  concept 
of  partial  Inhibition  of  convection,  shows  that  the  inhibition  mechanism 
can  yield  acceptable  distributions  of  surface  temperature.  The  results  of 
this  model  also  show  that:  (1)  the  edge  of  the  umbra  is  sharp,  even  for 
deep  spots;  (11)  deep  spots  produce  weak  bright  rings,  but  shallow  spots 
produce  Intense  bright  rings  in  conflict  with  observations;  and  (ill)  only 
a  shallow  surface  layer  of  the  sunspot  is  cool,  the  rest  being  warmer  than 
the  surroundings. 

Unbral  oscillations  in  sunspots  are  studied  and  identified  as  a 
resonant  response  of  the  umbral  atmosphere  to  forcing  by  oscillatory 
convection  in  the  subphotosphere.  The  full  linearized  equations  for 
magneto-atmospheric  waves  are  solved  nunerically  for  a  detailed  model  of 
the  umbral  atmosphere.  ^Resonant  "fast"  modes  are  found,  the  lowest  mode 
having  a  period  of  153  s^  typical  of  umbral  oscillations.  A  small  amount 
of  wave  energy  leaks  inti  the  corona  in  the  form  of  an  acoustic  wave  along 
the  jnafnet4oHLleld  litre . 

^ It  is  suggested  that  the  Sun's  radius  and  surface  temperature  vary 
with  the  solar  cycle  due  to  the  variation  of  total  magnetic  buoyancy  in 
the  convection  zone  over  the  cycle  of  the  solar  dynamo .  Thi s  mechanism 
might  produce  the  drop  in  surface  temperature  with  increasing  solar  activity 
observed  by  Livingston  (1978),  with  little  or  no  changf  n»  luminosity. 
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SUMMARY 


The  research  under  this  contract  has  focused  upon  dynamical  and  thermal 
phenomena  in  sunspots.  The  research  may  be  roughly  divided  into  three 
categories:  (i)  a  study  of  the  possible  mechanisms  for  cooling  a  sunspot, 
including  wave  cooling  and  inhibition  of  convection;  (ii)  a  study  of  resonant 
modes  of  magneto-atmospheric  waves  in  a  sunspot  umbra,  as  an  explanation  for 
umbral  oscillations;  and  (iii)  a  suggestion  that  the  Sun's  radius  and  surface 
temperature  may  vary  over  the  solar  cycle  due  to  changes  in  total  magnetic 
buoyancy  in  the  convection  zone.  The  last  topic,  while  somewhat  outside  of 
the  scope  of  the  initial  proposal,  nevertheless  arose  in  the  course  of  this 
research. 

The  question  of  the  mechanism  that  causes  a  sunspot  to  be  cooler  than 
its  surroundings  has  proved  to  be  quite  controversial.  The  traditional 
explanation  has  been  based  on  the  inhibition  of  convective  heat  transport  by 
the  strong  magnetic  field  in  the  sunspot  (Biermann  1941).  An  alternative 
mechanism,  based  on  cooling  by  an  outward  flux  of  Alfv€n  waves  (or  other 
magnetohydrodynami c  waves),  has  been  advocated,  most  recently  by  Parker 
(1974,  1975).  Both  of  these  mechanisms  have  been  investigated  further  under 
this  contract,  with  the  conclusion  that  the  inhibition  mechanism  is  by  far 
the  more  likely  explanation. 

Thomas  (1978;  section  2  of  this  report)  investigated  the  reflection  of 

Alfv6n  waves  in  the  umbral  atmosphere.  He  considered  the  reflection  of  waves 

propagating  upward  and  downward  in  the  umbra  from  a  generating  source 

(overstable  convection)  in  a  shallow  layer  just  below  the  umbral  photosphere. 

The  study  used  a  three-layer  model  of  the  umbral  atmosphere  that  reproduces 
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all  of  the  essential  features  of  the  vertical  structure  in  the  convection 
zone,  photosphere, chromosphere,  transition  region,  and  corona.  The  results 
show  very  strong  downward  reflection  of  A1 fvdn  waves  propagating  upward  into 
the  photosphere  and  low  chromosphere.  This  result  is  in  good  agreement  with 
the  observations  of  Beckers  (1976)  and  Beckers  and  Schneeberger  (1977),  which 
show  a  large  drop  in  wave  energy  density  from  the  photosphere  to  the  high 
chromosphere.  The  observations  put  an  upper  limit  on  the  upward  wave  energy 
flux  that  is  several  orders  of  magnitude  below  the  "missing"  flux  of  a 
sunspot.  Thus,  both  theory  and  observation  here  tend  to  rule  out  substantial 
cooling  by  upward-propagating  AlfvGn  waves. 

On  the  other  hand,  Thomas  (1978)  found  only  very  weak  upward  reflection 
of  Alfvdn  waves  propagating  downward  into  the  convection  zone.  This  left 
open  the  possibility  of  cooling  by  downward  propagating  waves.  Several 
difficulties  associated  with  cooling  by  downward -pro pa gating  waves  were 
noted,  however.  In  particular,  it  was  speculated  that  if  one  were  to  study 
more  realistic  wave  motions,  including  the  effects  of  compression  and 
buoyancy  (rather  than  the  pure  Alfv6n  mode),  then  one  would  find  much 
stronger  upward  reflection  in  the  convection  zone  due  to  the  rising  temperature 
(and  sound  speed)  with  depth.  Further  work  by  Scheuer  and  Thomas  (1980; 
section  5  of  this  report)  confirmed  this  expectation.  Using  the  same  model 
umbral  atmosphere,  but  considering  the  full  magneto-atmospheric  wave  equations 
(with  nonzero  horizontal  wavenumber),  Scheuer  and  Thomas  found  strong  upward 
reflection  in  the  convection  zone,  due  primarily  to  the  increasing  sound 
speed  with  depth.  This  strong  reflection  holds  true  even  for  waves  with 
horizontal  wavelengths  comparable  to  the  spot  diameter.  Thus,  cooling  by 
downward-propagating  waves  seems  to  be  ruled  out  too.  Our  conclusion  is 
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that  the  wave  cooling  mechanism  may  be  ruled  out  on  theoretical  and 
observational  grounds. 

In  order  to  assess  the  possibility  of  sunspot  cooling  by  the  mechanism 
of  inhibition  of  convection,  Clark  (1979,  Section  3  of  this  report)  studied 
simple  thermal  models  of  sunspots  which  involve  solutions  to  the  steady  heat 
equation  with  an  appropriate  distribution  of  thermal  conductivity.  Similar 
studies  had  been  made  by  Parker  (1974),  Eschrich  and  Krause  (1977),  and 
Spruit  (1977).  Clark  chose  a  simple  model  for  the  sunspot  which  allowed  him 
to  investigate  the  dependence  of  the  solution  on  the  various  spot  parameters, 
such  as  the  depth  of  region  of  inhibition.  He  found  that  his  model  sunspot, 
like  those  of  Eschrich  and  Krause  and  of  Spruit,  produced  surface  temperature 
distributions  very  much  like  that  observed  in  sunspots.  Detailed  results 
from  Clark's  model  include  the  following:  (i)  the  edge  of  the  umbra  is  sharp, 
even  for  deep  spots  (depth  of  inhibition  region  >  diameter);  (ii)  deep  spots 
produce  weak  bright  rings,  whereas  shallow  spots  produce  intense  bright  rings 
(which  aren't  observed);  and  (iii)  only  a  thin  surface  layer,  about  one 
temperature  scale  height  in  thickness,  is  cool,  with  the  deeper  parts  of  the 
spot  somewhat  hotter  than  the  surroundings.  The  three  thermal  models  of 
sunspots  due  to  Eschrich  and  Krause  (1977),  Spruit  (1977),  and  Clark  (1979) 
each  use  a  different  specification  for  the  thermal  conductivity  and  for  the 
spot  geometry,  and  yet  all  three  produce  acceptable  surface  temperature 
distributions.  This  lends  considerable  support  for  the  inhibition  mechanism. 

In  further  work  on  thermal  models  of  sunspots,  Clark  has  studied  the 
effects  of  a  depth-dependent  thermal  conductivity,  first  discussed  by  Spruit 
(1977).  Using  some  simple  analytical  solutions,  Clark  has  shown  that  the 
horizontal  spreading  of  the  heat  flux  depends  on  both  the  variable  conductivity 


7 


and  the  boundary  condition  at  the  bottom  of  the  convection  zone.  As  a 
particular  example,  consider  a  steady  heat  source  at  depth  zQ  below  the 
top  of  the  convection  zone,  and  take  the  convection  zone  to  be  infinitely 
deep.  At  the  top  of  the  convection  zone  assume  a  radiative  boundary  condition. 
If  the  thermal  conductivity  is  constant,  all  of  the  heat  emitted  by  the  source 
flows  out  the  top,  as  one  would  expect.  If,  however,  the  thermal  conductivity 
K  varies  with  depth  z,  the  result  can  be  very  different.  For  the 
conductivity  law  suggested  by  Spruit,  namely  K«(l+az)  with  a  =127  km, 
much  of  the  heat  emitted  by  the  source  flows  off  to  infinity  rather  than 
leaving  through  the  top.  In  fact,  the  fraction  of  the  heat  flux  leaving 
through  the  top  is 
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where  h  is  the  superadiabatic  temperature  scale  height  at  the  surface.  Thus 
for  deep  disturbances  (azQ»l)  very  little  of  the  released  heat  ever  reaches 
the  surface.  This  result  shows  that  the  boundary  condition  at  the  bottom  of 
the  convection  zone  is  likely  to  be  very  important. 

Scheuer  and  Thomas  (1980;  Section  5  of  this  report)  have  presented  a 

detailed  theory  of  umbral  oscillations,  identifying  them  as  the  lowest  resonant 

mode  of  fast  magneto-atmospheric  wave  in  the  umbral  atmosphere,  excited  by 

overstable  convection  in  the  subphotosphere.  Their  theory  is  similar  to  that 

of  Uchida  and  Sakurai  (1975),  except  that  Scheuer  and  Thomas  use  the  full 

magneto-atmospheric  wave  equations  rather  than  the  "quasi -Alfvdn"  approximation 

of  Uchida  and  Sakurai.  Numerical  solutions  of  the  full  wave  equations  are 

obtained  for  forced  and  free  oscillations  in  a  three-layer  model  of  a  sunspot 

umbra.  Several  new  features  of  the  wave  motion  arise  from  this  study,  leading 
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to  a  new  understandg  of  umbral  oscillations.  The  resonant  mode  is  trapped 
primarily  by  the  increasing  Alfvdn  speed  upward  into  the  chromosphere  and  by 
the  increasing  sound  speed  downward  into  the  convection  zone.  The  downward 
reflection  is  not  complete,  however;  a  small  fraction  of  the  total  energy 
escapes  to  large  heights  by  converting  into  the  form  of  a  pure  acoustic  wave 
along  the  vertical  magnetic  field  lines.  The  change  in  character  of  the  wave 
with  height,  from  Alfvdn-like  to  acoustic-like,  is  one  of  the  more  interest¬ 
ing  features  predicted  by  the  full  equations.  It  is  also  found  that  the 
increasing  sound  speed  in  the  convection  zone,  not  the  increasing  density 
(as  proposed  by  Uchida  and  Sakurai),  provides  the  necessary  upward  reflection 
for  the  trapped  mode.  For  numerical  values  of  the  model  parameters  based  on 
a  typical  sunspot,  the  lowest  resonant  fast  mode  has  a  period  of  153  s, 
typical  of  umbral  oscillations. 

On  another  topic,  Thomas  (1979,  Section  4  of  this  report)  has  suggested 

that  the  solar  radius  and  surface  temperature  may  vary  with  the  solar  cycle, 

due  to  changes  in  the  total  magnetic  buoyancy  in  the  convection  zone.  This 

idea  was  prompted  by  the  observations  of  Livingston  (1978),  which  show  a 

decrease  in  surface  temperature  of  the  Sun  with  increasing  solar  activity. 

Livingston  interprets  this  to  imply  a  corresponding  decrease  in  luminosity, 

but  Thomas  points  out  that  it  could  also  be  due  to  a  small  increase  in  the 

solar  radius  with  little  or  no  change  in  luminosity.  The  mechanism  proposed 

for  the  expansion  with  increasing  solar  activity  is  based  on  increasing 

magnetic  buoyancy  in  the  convection  zone  as  the  solar  dynamo  moves  toward 

solar  maximum.  Rough  estimates  of  the  difference  in  magnetic  flux  in  the 

convection  zone  between  solar  maximum  and  minimum  indicate  that  this  mechanism 

may  be  important.  The  historical  record  of  observations  of  the  solar  radius 
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is  quite  inconsistent,  but  there  is  evidence  of  variations  of  the  solar  radius 
with  the  solar  cycle,  with  maximum  radius  at  maximum  activity.  Accurate 
monitoring  of  the  solar  radius  over  a  solar  cycle  is  needed  to  evaluate  this 
idea  and  to  provide  needed  information  for  the  relation  between  solar 
luminosity,  temperature,  and  size. 
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ABSTRACT 

As  one  means  of  evaluating  the  possibility  that  sunspots  are  cooled  by  a  flux  of  Alfv6n  waves, 
the  reflection  of  vertically  propagating  AIfv6n  waves  in  a  three-layer  model  of  a  sunspot  umbra 
is  studied.  The  results  show  strong  downward  reflection  of  Alfv6n  waves  in  the  photosphere  and 
low  chromosphere,  with  very  little  wave  energy  penetrating  as  high  as  the  corona.  This  is  in 
agreement  with  recent  observations.  The  model  umbra  also  shows  very  weak  upward  reflection 
of  Alfv6n  waves  propagating  downward  into  the  convection  zone.  The  results  suggest  that,  if 
sunspots  are  indeed  cooled  by  Alfvdn  waves,  these  waves  must  escape  downward  into  the  solar 
interior. 

Subject  headings:  hydromagnetics  —  Sun:  sunspots 


I.  INTRODUCTION 

The  mechanism  that  causes  a  sunspot  to  be  cooler 
than  its  surroundings  is  still  not  understood.  Bier- 
mann’s  (1941)  suggestion  that  the  cause  is  the  in¬ 
hibition  of  convection  by  the  intense  magnetic  field 
in  the  sunspot  has  often  been  accepted  but  has  never 
developed  into  a  full-fledged  theory.  Parker  (1974a,  b , 
1975a,  b)  recently  revived  an  alternative  suggestion, 
based  on  earlier  work  by  others  (Danielson  1965; 
Musman  1967;  Savage  1969;  Moore  1973),  that  the 
cooling  is  due  to  a  flux  of  Alfv6n  waves  (or  other 
magneto-atmospheric  waves)  generated  by  overstable 
convection  in  the  umbral  subphotosphere  (see  also 
Roberts  1976).  It  is  possible,  of  course,  that  both 
mechanisms  operate  simultaneously;  indeed,  it  is  the 
inhibiting  effect  of  the  magnetic  field  on  convective 
motions  that  leads  to  the  overstability  that  supposedly 
generates  the  Alfv6n  waves. 

Parker’s  ideas  have  led  to  a  lively  controversy  over 
several  points.  Cowling  (1976)  has  objected  to  the 
high  thermal  efficiency  required  for  the  “refrigera¬ 
tion”  of  a  sunspot  by  a  flux  of  waves  (see  the  rebuttal 
by  Parker  1977).  By  means  of  a  simple  thermal  model, 
Parker  (1974a)  argued  that  the  inhibition  mechanism 
would  not  produce  the  observed  surface  temperature 
distribution  of  a  sunspot.  But  other  thermal  models 
(Eschrich  and  Krause  1977;  Spruit  1977)  have  pro¬ 
duced  acceptable  sunspots,  and  Clark  (1978)  has 
shown  that  this  point  hinges  on  what  one  assumes  the 
depth  of  a  sunspot  to  be  (see  also  Wilson  1971). 
Boruta  (1977)  has  argued  that  the  formation  of  a 
sunspot  is  more  readily  explained  by  the  wave-cooling 
mechanism.  Another  question  concerns  the  intrinsic 
stability  of  a  fully  developed  sunspot  (Parker  19756; 
Meyer,  Schmidt,  and  Weiss  1977),  which  may  depend 
on  the  cooling  mechanism. 

In  the  present  paper  we  deal  with  a  specific  problem 
associated  with  the  possible  cooling  by  A!fv6n  waves: 
the  escape  of  Alfv^n  waves  may  be  limited  by  re¬ 
flections  in  the  strong  vertical  density  gradients  in  the 


sunspot.  The  observations  of  Beckers  (1976)  and 
Beckers  and  Schneeberger  (1977),  taken  together,  in¬ 
dicate  a  strong  downward  reflection  of  Alfv6n  waves 
in  the  sunspot  atmosphere.  Beckers  (1976)  puts  an 
upper  bound  of  2.5  x  1010  ergs  cm-2  s"1  on  the 
upward  energy  flux  of  AIfv6n  waves  in  the  umbral 
photosphere,  and  Beckers  and  Schneeberger  (1977) 
put  an  upper  bound  of  4  x  107  ergs  cm- 2  s-1  on  the 
upward  energy  flux  of  Alfv6n  waves  in  the  corona 
above  a  sunspot  umbra.  The  latter  flux  limit  is  well 
below  the  “missing”  flux  of  5  x  1010  ergs  cm"2  s'1 
needed  to  cool  a  sunspot.  These  observations  leave 
open  the  possibility  that  the  Alfv6n  waves  escape 
downward  into  the  solar  interior. 

If  we  assume  a  purely  vertical,  uniform  umbral 
magnetic  field  B0  =  (0, 0,  B0),  in  a  Cartesian  co¬ 
ordinate  system  (x,  y ,  z)  with  z  upward,  and  assume 
purely  horizontal  motions  with  velocity  m  = 
[m(z,  /),  0, 0],  then  we  obtain  the  equation 


dt* 


=  vw  D 


<Pu 


(i) 


describing  vertically  propagating  pure  Alfv6n  waves 
in  the  umbra,  with  AIfv6n  speed 

Va(z)  -  [B02l4rrnp(z)]m 


that  varies  due  to  the  varying  density  p(z).  If  we 
assume  an  oscillatory  solution  of  the  form  u(z,  f)  = 
u(z)  exp  (iW ),  equation  (I)  becomes 


dH 
dz 2 


4* 


u  =  0. 


(2) 


For  a  specified  density  distribution  p(z)y  we  can  solve 
equation  (2)  to  determine  effective  reflection  coeffi¬ 
cients  for  Alfv6n  waves  in  the  atmosphere.  This  type 
of  problem  is  important  in  the  study  of  waves  in 
layered  media  (see  Brekhovskikh  1960),  and  is 
analogous,  for  example,  to  the  problem  of  reflection 
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of  light  in  a  medium  in  which  the  index  of  refraction 
varies  in  the  direction  of  propagation. 

Geronicolas  (1977)  has  considered  this  problem 
with  the  density  distribution1 

piz)  =  Poll  -  (Pi/Po)  tanh  kz]  (3) 

representing  the  sunspot,  where  p0»  Pi>  and  k  are 
parameters.  He  chose  values  of  the  parameters  to 
match  p(z)  with  Spruit’s  (1974)  model  of  the  convec¬ 
tion  zone.  The  level  z  =  0  corresponds  to  the  base  of 
the  umbra!  photosphere.  Geronicolas  found  only  a 
slight  net  downward  reflection  of  Alfv6n  waves  in  this 
atmosphere.  (To  be  precise,  he  found  that  if  Alfven 
waves  are  emitted  upward  and  downward  with  equal 
intensity  at  z  —  0,  then  the  net  reflection  is  such  that 
the  upward  energy  flux  is  46%  f  the  total,  instead  of 
50%  as  in  the  case  of  no  net  reflection.)  If  we  take 
this  result  together  with  the  observations  discussed 
above,  we  would  conclude  that  Alfven  waves  are  not 
cooling  sunspots. 

Objections  can  be  made  to  Geronicolas’s  density 
distribution  (3),  however.  It  gives  a  good  fit  to  den¬ 
sities  in  detailed  umbral  models  only  over  a  very 
limited  range  of  height  in  and  near  the  photosphere. 
The  density  (3)  approaches  uniform  values  a  short 
distance  above  and  below  the  photosphere,  whereas 
the  true  solar  density  continues  to  vary  strongly  at 
all  levels  above  and  below  the  surface.  The  net  effect 
of  this  inconsistency  is  that  Geronicolas's  model 
severely  underestimates  the  downward  reflection  of 
Alfven  waves  in  a  sunspot.  Here  we  propose  a  more 
realistic,  three-layer  model  of  the  vertical  density 
distribution  in  a  sunspot  umbra,  and  show  that  it 
leads  to  strong  downward  reflection  of  Alfven  waves. 

II.  THE  MODEL  UMBRA 

The  calculation  of  reflection  coefficients  for  waves 
in  a  variable  atmosphere  is  a  case  where  approximate 
methods,  such  as  the  WKB  method,  usually  fail.  For 
this  reason  it  is  desirable  to  have  a  model  umbra  that 
is  simple  enough  to  allow  exact  solution  of  the  wave 
equation  (2),  but  detailed  enough  to  give  a  fairly 
accurate  representation  of  the  vertical  sunspot 
structure. 

We  adopt  as  our  model  umbra  an  atmosphere  with 
a  uniform  vertical  magnetic  field  and  the  three-layer 
temperature  distribution  shown  in  Figure  1.  Layer  2 
is  isothermal  at  temperature  T0  and  represents  the 
umbral  photosphere  and  chromosphere.  The  density 
in  layer  2  is  thus  given  by 

p(z)  =  p0  exp  ( -zjHc)  f  0  <  z  <  2t ,  (4) 

where  H0  -  RT0lg  is  the  density  scale  height  in  this 
layer.  Layer  3,  which  represents  the  corona,  is  iso¬ 
thermal  at  temperature  T,  ( >  T0),  and  the  transition 

1  The  formula  for  #><r)  has  been  changed  here  to  agree  with 
our  present  sign  convention;  Geronicolas  took  the  ?-axis  to 
point  downward.  The  problem  of  Alfven  wave  reflection  with 
the  iame  density  distribution,  (3),  has  been  considered  by 
Adam  (1976)  in  a  different  context 


Tf  z  > 


Fig.  1. — Schematic  diagram  of  the  three-layer  temperature 
distribution  in  the  model  umbra. 


region  is  represented  as  a  discontinuity  in  tempera¬ 
ture  (and  density)  at  z  -  zt.  Pressure  must  be  con¬ 
tinuous  across  z  =  zt,  however,  and  thus  the  density 
in  layer  3  is  given  by 

p(z)  =  Po  exp  (-z,///0)^  exp  [-(z  -  z,)///J , 

z  >  z, ,  (5) 

where  Hv  =  RTJg  is  the  scale  height  in  this  layer. 

Layer  I,  which  represents  the  umbral  convection 
zone,  is  assumed  to  have  a  linear  temperature  distribu¬ 
tion  of  the  form  T(z)  =  T0  —  P z,  with  P  >  0.  The 
value  of  p  may  be  chosen  equal  to  or  slightly  greater 
than  the  adiabatic  gradient  p,  =  g/cp.  The  corre¬ 
sponding  density  distribution  in  layer  1  is  then 

\  (l  —  o)la 

1  -  ,  z  <  0  ,  (6) 

where  a  =  Rpjg  is  a  dimensionless  measure  of  the 
temperature  gradient  p . 

Layers  2  and  3  are  the  same  as  the  two-layer  model 
umbra  studied  by  Uchida  and  Sakurai  (1975)  in  con¬ 
nection  with  umbral  oscillations.  The  same  tempera¬ 
ture  structure  as  in  layers  2  and  3,  but  with  a  vertical 
dipole  magnetic  field,  was  studied  by  Uchida  and 
Kaburaki  (1974)  in  connection  with  the  heating  of 
the  chromosphere  and  corona  above  sunspots.  We 
shall  see  here  that  layer  3  (the  corona)  is  unimportant 
as  far  as  the  reflection  of  Alfv6n  waves  that  might  be 
cooling  a  sunspot  is  concerned. 

III.  DOWNWARD  REFLECTION  OF  ALFVEN  WAVES 

We  imagine  that  Alfven  waves  are  generated  in  a 
thin  layer  near  z  =  0  in  the  model  umbra,  with  equal 
intensity  upward  and  downward,  and  we  then  ask 
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how  much  of  the  wave  energy  is  reflected  downward 
or  upward.  First  we  consider  the  downward  reflection 
from  layers  2  and  3.  To  isolate  this  aspect  of  the 
problem,  we  imagine  for  the  moment  that  the  region 
z  <  0  is  homogeneous  with  density  p0  and  that  there 
is  a  train  of  plane  Alfv^n  waves  of  fixed  frequency 
and  unit  amplitude  incident  upon  z  =  0  from  below. 
We  then  write,  for  z  <  0, 

(7) 

where  v*-  =  (2?03/4ffMPo)1/a  is  the  Alfv6n  speed  at 
z  =  0.  The  first  term  on  the  right-hand  side  of 
equation  (7)  is  an  incident  (upward  propagating)  wave 
of  unit  amplitude.  The  second  term  on  the  right-hand 
side  is  the  reflected  (downward  propagating)  wave  of 
undetermined  amplitude  U .  The  reflected  amplitude  is 
determined  by  matching  the  solution  (7)  to  the  ap¬ 
propriate  solution  in  layers  2  and  3.  The  proper 
matching  conditions  at  the  interfaces  (z  =  0  and 
z  =  zt)  are  that  the  velocity  u  and  the  x-component  of 
the  magnetic  field  perturbation  be  continuous  across  the 
interface;  the  latter  condition  is  equivalent  to  requiring 
dujdz  to  be  continuous  across  the  interface.  Once  U 
is  determined,  the  downward  reflection  coefficient  for 
energy,  Kd,  is  given  by  Rd  =  \U\2. 

First,  let  us  ignore  layer  3  and  assume  that  the 
density  distribution  of  layer  2  (eq.  [4])  extends  to 
z  =  +oo.  (The  reason  for  doing  this  will  become  clear 
in  a  moment.)  As  first  shown  by  Ferraro  (1954;  see 
also  Ferraro  and  Plumpton  1958),  the  wave  equation 
(2)  with  density  distribution  (4)  can  be  transformed 
into  Bessel’s  equation. 


r2d2u  rdu  t  _ 

c V  +  C5j  +  C“  =  0> 


(8) 


with  the  change  of  independent  variable  i  = 
a  exp  (-z/2//0),  where  a  -  2H0w/vAo.  Note  that  the 
parameter  a  can  also  be  written  as  a  ~  4tt//0/A,  where 
A  =  2nvAJco  is  the  wavelength  of  the  incident  Alfv6n 
wave;  thus  a  is  a  nondimensional  wavenumber.  The 
general  solution  of  equation  (8)  is  written  as  usual  as 


m  =  AMO  +  BY0(O.  (9) 

We  require  the  solution  to  be  finite  as  z  -*  +  oo  (i.e,, 
outgoing  wave  as  z  ^  +  oo),  so  B  =  0.  Then,  matching 
u  and  dtydz  in  solutions  (7)  and  (9)  at  z  —  0  ({  —  a), 
we  obtain 


and 


U<*)  +  Ux{a)' 


V  _  \Jxia)  ±  jUa)  1 . 

Vi(fl)  -  ^o(fl)J 


(10) 

(11) 


Thus,  Rd  =  |l/|a  =  I  for  any  value  of  a,  and  there  is 
total  downward  reflection  of  Alfv6n  waves  of  any 


wavelength.  This  result  should  not  be  surprising;  any 
density  distribution  with  p-*0asz-*  +  oo  will  give 
/**  -  1.  Of  more  importance  here  is  the  fact  that,  for 
wavelengths  expected  for  Alfv6n  waves  associated 
with  sunspot  cooling,  the  reflection  occurs  quite  low 
in  the  umbral  atmosphere.  The  Alfv6n  waves  do  not 
propagate  far  enough  upward  to  cool  the  spot. 

A  reasonable  estimate  of  the  scale  height  in  layer  2 
is  H0  =  100  km  corresponding  roughly  to  T0  =*= 
4500  K.  Beckers  (1976)  estimates  that  upward-prop¬ 
agating  Alfv6n  waves  that  might  be  cooling  a  sunspot 
would  have  wavelengths  in  the  range  15-40  km.  For 
these  wavelengths,  a  £  10ir,  and  we  can  use  the 
asymptotic  forms 

■/»(a)  =  (^)1,acos(fl-S)’ 

■/i(a)  =  (^)1'2sin  (a-5)- 

Then  equations  (9)  and  (10)  give 

1*1  «  (2ira)1/3y0({) . 

As  z  +co({  -►  0),  y0({)  I.  and  the  wave  amplitude 
approaches  a  uniform  value  (2 wa)1/a,  while  the  kinetic 
energy  density  decreases  rapidly  because  of  the 
exponentially  decreasing  mass  density.  (The  magnetic 
energy  density  also  decreases  exponentially  with  in¬ 
creasing  z  [Ferraro  1954).)  We  can  write 

p|«|2  -  2napJ02(0  <  2napQ  exp  ( —  z///0) , 

and  thus  most  of  the  wave  energy  is  reflected  down¬ 
ward  within  the  first  few  scale  heights  above  z  =  0. 
A  reasonable  estimate  of  the  height  of  the  transition 
region  is  zf  =  2000  km.  Since  this  is  20  scale  heights 
above  z  =  0  in  our  model,  only  a  very  small  fraction 
of  the  wave  energy  penetrates  up  to  the  transition 
region,  and  the  effect  of  the  corona  (layer  3)  is  un¬ 
important  as  far  as  the  reflection  of  the  Alfv^n  waves 
is  concerned.  For  completeness  the  solution  with 
layer  3  included  is  presented  in  the  Appendix.  The 
downward  reflection  coefficient  is  again  unity,  and 
the  results  of  the  present  section  are  shown  to  be  valid 
in  the  appropriate  limit. 

Note  that  the  solution  d({)  oc  J0(0  oscillates  with 
increasing  z  for  z  >  0  ({  <  a)  before  becoming  essen¬ 
tially  uniform  as  z  -*  +ao  (£  0).  Since,  for  example, 

yo(0.4)  =  0.96,  the  wave  amplitude  reaches  %70  of 
its  asymptotic  value  at  height  z0  =  -2 H0  log  (0.4/a). 
The  height  z0  is  greater  for  waves  of  shorter  wave¬ 
length  (larger  a);  waves  of  shorter  wavelength  (higher 
frequency)  propagate  farther  in  the  direction  of  in¬ 
creasing  Alfv6n  speed.  But  even  for  a  wavelength  as 
short  as  10  km  (a  =  40w),  the  height  z0  =  1150  km 
is  well  below  the  transition  zone.  For  z  >  z0,  the  wave 
amplitude  is  essentially  constant,  and  the  energy 
density  continues  to  decrease  exponentially  in  pro¬ 
portion  to  the  mass  density.  The  result  that  the  strong 
downward  reflection  of  Alfv6n  waves  occurs  very  low 
in  the  sunspot  atmosphere  helps  to  justify  our  initial 
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assumption  of  a  purely  vertical  magnetic  field  in  the 
umbra. 

In  the  observations  of  Beckers  (1976)  and  Beckers 
and  Schneeberger  (1977),  the  directly  observed  quan¬ 
tity  is  the  rms  velocity  of  Alfv6n  waves  causing  the 
nonthermal  line  broadening.  The  values  they  obtain 
for  the  rms  velocity  are  roughly  1.5  km  s"1  for  the 
photosphere  and  7.5  km  s"1  for  the  corona.  This 
small  observed  change  in  rms  velocity  with  height  is 
consistent  with  the  solution  presented  here. 

IV.  UPWARD  REFLECTION  OF  ALFV^N  WAVES 

To  isolate  the  problem  of  determining  the  upward 
reflection  coefficient  in  layer  1,  we  imagine  the  region 
z  >  0  to  be  homogeneous  with  density  p0  and  con¬ 
sider  a  train  of  plane  Alfv^n  waves  of  fixed  frequency 
incident  upon  2  =  0  from  above.  We  write,  for  z  >  0, 

0 -«*[*.(. 

(12) 

The  first  term  on  the  right-hand  side  of  equation  (12) 
is  an  incident  (downward  propagating)  wave  of  unit 
amplitude,  and  the  second  term  on  the  right-hand 
side  is  the  reflected  (upward  propagating)  wave  of 
undetermined  amplitude  K 

The  basic  wave  equation  (2)  can  be  solved  exactly 
in  layer  1,  where  the  density  is  given  by  equation  (6). 
If  we  let 


V  = 


a 

a  +  I 


02 

mJ 


where  c  »  2//0oj/yAo  as  before,  then  the  general  solu¬ 
tion  of  equation  (2)  in  layer  I  can  be  expressed  as 
(Watson  1966) 


u(v)  =  VM// “>( T,)  +  BHv'*\ri)] ,  (13) 

where  v  =  a/(a  +  1)  and  where  //VU)  and  Hv{2)  are 
Hankel  functions  of  the  first  and  second  kinds  of 
order  v.  For  z  -oo  (17  +  00),  the  principal  asymp¬ 

totic  forms  of  the  Hankel  functions  are 


(i + »)-]}• 

H-“w ~  -  (l + j)*]}  • 

Thus,  as  z  ->  -00  (r)  ->  +  00),  we  see  that  repre¬ 
sents  a  wave  in  the  +2  (-17)  direction  and  Hv{2) 
represents  a  wave  in  the  -z  (  +  »j)  direction.  Here  we 
want  only  an  outgoing  wave  as  2  -*  -00,  so  we  must 
take  A  =  0  in  equation  (13). 

If  we  now  match  ft  and  dft/dz  in  solutions  (12)  and 
(13)  at  2  =  0,  we  obtain 


1*  //  (»>/-  \  _l_  iff  13)/  0  \ 

/U“+iJ 

>r’  v°  +  i/  ^  \ZTi)\ 

04) 


and 

M.-fi)  -  "•-■“tfi)]-  ,l5» 

The  Hankel  functions  are  complex  valued,  so  )  V\  ^  l 
here.  To  evaluate  V  we  must  first  choose  an  appro¬ 
priate  value  of  a,  the  dimensionless  temperature 
gradient  in  layer  1.  To  avoid  tedious  computation  we 
shall  choose  values  of  a  such  that  v  =  a /(a  -1-  1)  and 
v  -  1  =  —  1  /(«  -h  1)  are  among  the  fractional  orders 
for  which  tables  of  Bessel  functions  are  available.  The 
values  a  =  \  (giving  v  -  v  -  I  =  -  J)  and  a  =  $ 
(giving  v  -  l  and  v  —  1  =  -$)  seem  most  appro¬ 
priate.  We  can  write 

a  =  —  =  (Cp  Cv\  Cp^  —  / V  1\  £ 
g  \  cp  )  g  \  y  /  ft  ’ 

where  ft  =  g/c„  is  the  adiabatic  temperature  gradient. 
Thus, 


(16) 


We  can  place  bounds  on  a  as  follows.  We  want 
0/ft  >  1,  which  requires  a  >  (y  -  1  )/y.  Most  models 
of  the  umbral  subphotosphere  show  density  increasing 
with  depth  more  rapidly  than  a  linear  variation,  and 
we  see  from  equation  (6)  that  this  requires  a  5 
Thus,  we  have  (y  —  l)/y  <<*<■£.  The  lower  bound 
equals  f  for  y  =  j  or  $  for  y  =  f .  Thus,  the  two  values 
a  =  i  and  a  =  $  pretty  well  cover  the  expected  range. 

Values  of  V  and  of  the  upward  reflection  coefficient 
(for  energy)  RM  =  \  V\ 2  have  been  obtained  by  ex¬ 
pressing  the  Hankel  functions  in  equation  (15)  in 
terms  of  Bessel  functions  of  the  first  kind,  and  using 
tabled  values  of  the  Bessel  functions  of  fractional 
order.  Table  1  gives  values  of  Ru  for  different  dimen¬ 
sionless  wavenumbers  a.  For  a  »  1,  corresponding  to 
the  short  wavelengths  (15-40  km)  discussed  in  the 
previous  section,  we  can  find  a  simple  asymptotic 
expression  for  Using  Hankers  asymptotic  ex¬ 
pansion 


X  exp  |  -fjr,  -  ^  +  ^*||  •  Oj  »  0 . 


we  find,  to  first  order  in  a'1, 


V  = 


i(l  -  «) 
4a 


R. 


IH* 


(i  -«) 
16? 


a 


a  1  . 


(17) 


and  hence 
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a  «  =  i  «  =  i 


«  +  1  XlHo  Ru  XIHq  Ru 


10 .  8.38  X  10’1  6.9  X  10*  9.42  x  10’1  1.6  x  10’4 

5 .  1.68  2.7  x  10“4  1.88  6.3  x  10'4 

1 .  8.38  5.0  x  10’*  9.42  1.1  x  10’1 

5  x  10'1 .  1.68  x  101  1.4  x  10"a  1.88  x  101  3.2  x  10"* 

I  x  10'1 .  8.38  x  10*  8.0  x  10’*  9.42  x  W  1.7  x  10’1 

5  x  10  * .  1.68  x  10*  1.4  x  10’x  1.88  x  10*  2.8  x  10’1 

1  x  10'* .  8.38  x  10*  3.1  x  10’1  9.42  x  10*  5.6  x  10’1 

5  x  10'3 .  1.68  x  103  4.0  x  10  l  1.88  x  10*  6.6  x  10* 

1  x  10'3  .  8.38  x  103  5.8  x  10’1  9.42  x  10*  8.3  x  10“l 


This  expression  is  accurate  to  two  significant  figures 
for  a  >  5. 

Table  1  and  expression  (17)  show  that  the  upward 
reflection  from  the  umbral  subphotosphere  is  very 
weak.  For  A  =  40  km  (a  —  lOrr)  and  a  =  equation 
(17)  gives  Ru  -  1.58  x  10  “5.  Of  course,  we  cannot 
put  any  observational  limit  on  the  wavelength  of 
Alfv6n  waves  that  might  be  propagating  downward 
into  the  Sun.  But  even  if  these  waves  had  wavelengths 
comparable  to  the  size  of  the  whole  umbra,  the 
reflection  coefficient  is  still  small  (Ru  =  0.175  for 
a  =  $,  A  =  94.25//0  =  9425  km). 

V.  DISCUSSION 

The  model  umbra  considered  here  predicts  total 
downward  reflection  and  negligible  upward  reflection 
of  Alfv6n  waves  in  a  sunspot.  The  downward  reflection 
in  the  photosphere  and  chromosphere  is  so  strong  that 
very  little  wave  energy  reaches  as  high  as  the  transition 
region  and  corona.  Thus,  if  Alfv&i  waves  arc  indeed 
cooling  a  sunspot,  we  conclude  that  these  waves  must 
propagate  downward  into  the  solar  interior.  This 
same  conclusion  was  reached  by  Wilson  (1975)  and 
by  Beckers  and  Schneeberger  (1977)  on  the  basis  of 
observations. 

The  very  small  upward  reflection  coefficient  for 
Alfv£n  waves  in  our  model  umbra  indicates  that 
cooling  by  downward-propagating  waves  is  possible. 
However,  because  of  great  uncertainties  about  the 
structure  of  a  sunspot  below  the  solar  surface,  this 
result  can  be  considered  only  as  suggestive.  If  the  spot 
is  fairly  deep  and  the  magnetic  flux  rope  nearly  vertical 


in  the  convection  zone,  as  suggested  by  Meyer  et  at . 
(1974),  then  the  present  model  is  fairly  accurate.  In 
order  for  the  Alfv6n  waves  to  cool  the  sunspot,  they 
must  propagate  sufficiently  far  from  the  surface  layers 
of  the  sunspot  before  their  energy  is  dissipated  into 
heat  and  thoroughly  mixed  to  become  part  of  the 
overall  radial  heat  flux  in  the  solar  interior.  This 
process  could  be  limited  by  scattering  and  dissipation 
of  the  Alfv6n  waves  due  to  inhomogeneities  in  the 
density  and  magnetic  field  in  the  subsurface  structure 
of  the  sunspot.  Much  more  work  needs  to  be  done  to 
fairly  evaluate  the  possibility  of  cooling  by  downward- 
propagating  Alfv6n  waves. 

It  should  also  be  noted  that  the  pure  Alfvta  waves 
considered  here  are  a  very  special  type  of  wave, 
involving  only  horizontal  motions.  Umbral  convection 
may  not  produce  this  kind  of  motion  very  efficiently, 
and  one  should  perhaps  consider  more  general  mag¬ 
neto-atmospheric  waves  with  vertical  motions  and 
compression  as  well.  But  these  more  general  waves 
will  no  doubt  suffer  much  greater  upward  reflection 
than  the  pure  Alfv6n  waves  because  of  the  increasing 
sound  speed  and  rapidly  changing  buoyancy  force  in 
the  umbral  subphotosphere.  This  presents  an  addi¬ 
tional  difficulty  for  the  wave-cooling  hypothesis. 

I  am  indebted  to  Jacques  Beckers,  Alfred  Clark, 
Jr.,  and  John  Molyneux  for  helpful  discussions.  Most 
of  this  work  was  done  during  a  visit  at  Sacramento 
Peak  Observatory.  This  work  was  supported  by  the 
Air  Force  Geophysics  Laboratory  under  contract 
F 19628-77-C-0079. 


APPENDIX 


Here  we  present  the  calculation  of  the  downward  reflection  coefficient  Rd  in  the  case  where  layer  3  (the  corona) 
of  the  model  umbra  is  included.  A  similar  analysis  has  been  given  by  Hollweg  (1972)  in  a  different  context.  The 
solution  of  the  wave  equation  (2)  in  layer  3  that  represents  an  outgoing  wave  as  z  +oo  is 


where 


*  =  C7o(f), 


(  =  frexp[-(z-z<)/2/7l]. 


«>A. 


=  m 


cxp(-Zt/2H0), 


(Al) 

(A2) 


280 


THOMAS 


and  Hi  =  RTJg  is  the  scale  height  in  layer  3.  We  have  solution  (7)  in  layer  1,  solution  (9)  in  layer  2  (with 
B  #  0),  and  solution  (Al)  in  layer  3.  Matching  u  and  dujdz  across  z  =  0  and  z  =  z(,  we  obtain  (after  some 
algebraic  manipulation) 


,  2[*Ji(b)YMt)  - 
M  +  iN 


(A3) 


and 


„  2[7,(b)/1({,)  -  sJdbVMM 

B  M  +  iN 

c  2[r,(C)/.({,)  -/,({,)  r,(oi 

M  +  iN 


U 


M  -iN 
M  +  iN’ 


(A4) 

(A5) 

(A6) 


where  =  a  exp  (-z,/2//0)  and 


*  =  {^o(h)[y»(W  Ko(fl)  -  Ki(tVo(«)]  -  sJiWUl,)  Ua)  -  r0(C,Vc(a)T} ,  (A7) 

N  =  {/„(*)[/,(£,)  r,(fl)  -  K, («/,(«)]  -  sJi(b)V0&)  7,(0)  -  Ko&Vito]}  •  (A8) 

Since  M  and  N  are  real,  we  have  R„  =  |t/|*  =  1  for  all  values  of  the  parameters  a,  b,  and 
Considerable  simplification  of  the  coefficients  A,  B,  and  C  is  possible  for  our  purposes.  If  we  adopt  the  values 
T0  =  4500  K,  Ti  =  10*  K,  H0  =  100  km,  z,  =  2000  km,  then  s  =  3  x  10'®,  and  it  is  sufficient  to  take  s -*■  0  in 
equations  (A3HA8).  The  parameter  is  also  small  as  long  as  a  is  not  too  large  ({,  <  5.7  x  10'*  for  A  >  10  km), 
and  we  can  also  take  {,  -»  0  in  equations  (A3MA8),  giving 

2  2 

A  =  Ua)  +  UM)  ’  5==0’  C  =  Ub)[J0(a)  +  i/Ja)}  *  (A9) 

In  this  limit  the  solution  in  layer  2  is  the  same  as  that  obtained  in  §  III  where  we  ignored  layer  3.  Also,  for 

A  ^  10  km,  b  <  0.085  and  J0(b)  »  1 .  Thus,  the  only  important  effect  of  including  layer  3  for  the  parameter  values 

assumed  here  is  to  change  the  length  scale  in  the  solution  above  z  =  zt. 


REFERENCES 


Adam,  J.  A.  1976,  J.  Phys.  Ay  9,  L193, 

Beckers,  J.  M.  1976,  Ap .  J 203,  739. 

Beckers,  J.  M.,  and  Schneeberger,  T.  J.  1977,  Ap.  7.,  215,  356. 
Biermann,  L.  1941,  Vierteljahrsschr .  Astr.  Geselisch 76,  194. 
Boruta,  N.  1977,  Ap.  7.,  215,  364. 

Brekhovskikh,  L.  M.  1960,  Waves  in  Layered  Media  (New 
York:  Academic  Press). 

Clark,  A.,  Jr.  1978,  in  preparation. 

Cowling,  T.  G.  1976,  M.N  R.A.S.,  177,  409. 

Danielson,  R.  E.  1965,  in  I AU  Symposium  No.  22 y  Stellar  and 
Solar  Magnetic  Fields ,  ed.  R.  Ltist  (Amsterdam:  North- 
Holland),  p.  315, 

Eschrich,  K.-O.,  and  Krause,  F.  1977,  Astr.  Nach.y  298,  1. 
Ferraro,  V.  C.  A.  1954,  Ap.  7.,  119,  393. 

Ferraro,  V,  C.  A.,  and  Plumpton,  C.  1958,  Ap.  7.,  127,  459. 
Geronicolas,  E.  A.  1977,  Ap.  7.,  211,  966. 

Hollweg,  J.  V.  1972,  Cosmic  Electrodyn.,  2,  423. 

Meyer,  F.,  Schmidt,  H.  U.,  and  Weiss,  N.  O.  1977, 
M.N.R.A.S. » 179,  741. 

Meyer,  F.,  Schmidt,  H.  U.,  Weiss,  N.  O.,  and  Wilson,  P.  R. 
1974,  M.N.R.A.S.,  169,  35. 


Moore,  R.  L.  1973,  Solar  Phys.9  30,  403. 

Musman,  S.  1967,  Ap.J.,  149,  201. 

Parker,  E.  N.  1974a,  Solar  Phys.,  36,  249. 

- .  19746,  Solar  Phys.,  37,  127. 

- .  1975a,  Solar  Phys.,  40,  275. 

- .  19756,  Solar  Phys.,  40,  291. 

- .  1977,  M.N.R.A.S.,  179,  93p. 

Roberts,  B.  1976,  Ap.  7.,  204,  268. 

Savage,  B.  D.  1969,  Ap.  7.,  156,  707. 

Spruit,  H.  C.  1974,  Solar  Phys.,  34,  277. 

- .  1977,  Solar  Phys.,  55,  3. 

Uchida,  Y.,  and  Kaburaki,  O.  1974,  Solar  Phys.,  35,  451. 
Uchida,  Y.,  and  Sakurai,  T.  1975,  Pub.  Astr.  Soc.  Japan ,  27, 
259. 

Watson,  G.  N.  1966,  A  Treatise  on  the  Theory  of  Bessel 
Functions  (2d  ed.;  Cambridge:  Cambridge  University 
Press),  p.  97. 

Wilson,  P.  R.  1971,  Solar  Phys.,  21,  101. 

- .  1975,  Solar  Phys.,  42,  333. 


John  H.  Thomas:  Department  of  Mechanical  and  Aerospace  Sciences,  University  of  Rochester,  Rochester, 
NY  14627 


THERMAL  MODELS  OF  SUNSPOTS 


ALFRED  CLARK,  JR. 

Department  of  Mechanical  and  Aerospace  Sciences  and  C.  E.  Kenneth  Mees  Observatory,  University  of 
Rochester ,  Rochester,  N.  Y.  14627 ,  U  S  A. 

(Received  20  July;  in  revised  form  3  December,  1978) 

Abstract.  We  study  simple  thermal  models  of  sunspots  based  on  the  concept  of  partial  inhibition  of 
convection  by  strong  magnetic  fields.  As  in  other  similar  studies  (Parker,  1974a;  Eschrich  and  Krause, 
1977;  Spruit,  1977a),  the  calculations  involve  solutions  of  the  heat  equation  with  an  appropriate 
distribution  of  thermal  conductivity. 

The  simplicity  of  the  present  model  allows  a  detailed  study  of  the  dependence  of  the  solution  on  the  spot 
parameters,  such  as  the  depth  of  the  region  in  which  convection  is  inhibited.  The  most  important  specific 
results  from  the  model  are:  (1)  the  edge  of  the  umbra  is  sharp,  even  for  deep  spots;  (2)  deep  spots  produce 
weak  bright  rings  in  the  surrounding  atmosphere,  whereas  shallow  spots  produce  intense  rings  which  are 
difficult  to  reconcile  with  observations;  (3)  only  a  surface  layer  of  a  spot,  with  thickness  of  the  order  of  the 
temperature  scale  height,  is  cool. 

The  present  model,  like  those  of  Eschrich  and  Krause  (1977)  and  Spruit  (1977a)  yields  surface 
temperature  distributions  resembling  sunspots.  Since  the  three  models  all  use  different  descriptions  of  the 
convective  heat  transport,  we  conclude  that  the  major  predictions  of  the  inhibition  theory  are  relatively 
insensitive  to  model  details. 

1.  Introduction 

Many  years  ago,  Biermann  (1941)  suggested  that  a  sunspot  is  dark  because  the  strong 
magnetic  field  of  the  spot  suppresses  the  convection  that  normally  transports  the 
solar  flux.  Discussions  by  Hoyle  (1949)  and  Cowling  (1953)  called  attention 
to  the  possibility  that  the  magnetic  field  reduces,  but  does  not  suppress  entirely,  the 
convective  heat  flux.  Later  magnetohydrostatic  models  of  sunspots  by  Chitre  (1963), 
Deinzer  (1965),  Chitre  and  Shaviv  (1967),  Yun  (1970),  Mullan  (1973),  and  Busse 
(1973)  incorporated  this  concept  in  various  ways. 

Interest  in  the  problem  has  been  revived  by  Parker’s  (1974a,  b,  1975a,  b)  vigorous 
advocacy  of  an  alternative  mechanism,  originally  suggested  by  Danielson  (1965), 
whereby  sunspots  are  cooled  by  the  radiation  of  Alfv6n  waves.  A  number  of  recent 
papers,  in  addition  to  Parker’s,  have  presented  calculations  bearing  on  the  AIfv£n 
wave  mechanism  (e.g.,  Roberts,  1976;  Boruta,  1977;  Thomas,  1978),  or  the 
inhibition  of  convection.  The  subject  is  controversial  (Cowling,  1976a,  Parker,  1977) 
and  most  authors  have  dealt  with  one  or  the  other  of  the  two  mechanisms  (the  present 
work  is  no  exception).  However,  it  is  worth  emphasizing  that  they  are  not  mutually 
exclusive.  It  may  well  be  that  both  are  important  in  sunspots.  The  work  reported  here 
is  concerned  entirely  with  the  mechanism  of  inhibition  of  convection,  so,  apart  from 
an  occasional  comment,  we  do  not  discuss  further  the  cooling  by  AIfv6n  waves. 

The  model  presented  here  is  related  to  the  models  of  Parker  (1974a),  Eschrich  and 
Krause  (1977),  and  Spruit  (1977a).  All  four  models  are  appropriately  called  thermal 
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models,  because  they  involve  the  calculation  of  a  temperature  distribution  from 
some  form  of  heat  transport  equation.  These  models  do  not  deal  with  the  mechanical 
aspects  of  a  sunspot.  Thus  they  provide  a  necessary  but  not  sufficient  test  for  a 
proposed  cooling  mechanism.  The  additional  constraints  imposed  by  mechanical 
equilibrium  may  rule  out  mechanisms  which  are  allowed  by  thermal  models. 
Nevertheless,  the  models  are  very  useful  in  studying  the  thermal  balance  of  a 
sunspot. 

The  four  thermal  models  under  present  discussion  are  all  based  on  the  assumption 
that  the  turbulent  convective  heat  transport  can  be  modeled  by  some  kind  of 
conductive  law.  There  is  no  unique,  or  even  generally  accepted  way  to  do  this,  and,  in 
fact,  such  an  assumption  is,  with  present  knowledge,  un verifiable.  Because  of  this, 
one  must  accept  at  the  outset  the  impossibility  of  absolutely  verifying  or  rejecting  the 
inhibition  concept  on  the  basis  of  thermal  models.  The  proper  objective  is  much 
more  modest:  namely,  an  attempt  to  assess  plausibility.  If  a  number  of  different  ways 
of  describing  the  turbulent  convective  heat  transport  lead  to  models  which  look  like 
sunspots,  then  the  magnetic  inhibition  concept  becomes  more  plausible.  If,  on  the 
other  hand,  only  very  special  circumstances  can  yield  a  model  resembling  a  spot,  then 
some  doubt  is  cast  on  the  role  of  magnetic  inhibition. 

Consider  now  in  more  detail  the  nature  of  the  four  thermal  models.  In  each  case, 
one  solves  the  steady-state  heat  equation  in  a  half-space,  with  the  heat  flux 
asymptotically  equal  to  the  solar  flux  deep  below  the  surface,  and  with  an  appropriate 
boundary  condition  (usually  the  Stefan -Boltzmann  law)  on  the  upper  surface.  The 
success  of  the  model  is  then  determined  by  comparing  the  calculated  distribution  of 
surface  temperature  with  temperature  distributions  observed  in  sunspots.  In  con¬ 
structing  a  model,  there  are  three  major  choices:  (i)  the  geometry  of  the  region  of 
strong  magnetic  field  where  the  convective  heat  transport  is  affected  (for  brevity,  we 
will  call  this  region  the  spot  in  what  follows);  (ii)  the  form  of  the  conductive  law 
outside  the  spot;  (iii)  the  alteration  of  the  conductive  law  within  the  spot.  We  now 
summarize  briefly  how  these  three  choices  were  made  in  the  four  models  under 
discussion. 

Eschrich  and  Krause  (1977)  chose  a  right  circular  cylinder  for  the  spot  shape. 
Outside  of  the  spot,  they  used  a  Fourier  law  for  the  heat  flux  with  an  isotropic, 
constant  conductivity.  Within  the  spot  they  used  a  constant  anisotropic  conductivity, 
with  the  vertical  conductivity  (along  the  magnetic  field)  being  considerably  less  than 
the  horizontal  conductivity. 

Spruit  (1977a)  also  used  a  cylindrical  spot  geometry.  He  used  a  more  complicated 
description  of  the  heat  transport,  based  on  mixing  length  theory.  His  conductivity 
was  both  depth -dependent  and  anisotropic.  Outside  the  spot,  the  anisotropy  was 
small.  Inside  the  spot,  Spruit  took  a  vanishing  horizontal  conductivity.  Thus  although 
Eschrich  and  Krause  (1977)  and  Spruit  (1977a)  agree  that  the  convective  transport  is 
anisotropic,  they  do  not  agree  on  the  sense  of  the  anisotropy.  As  we  discussed  earlier, 
this  kind  of  ambiguity  is  inherent  in  any  attempt  to  model  turbulent  convective 
transport. 
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The  present  work  fits  into  the  catalog  of  models  between  the  two  works  discussed 
above,  since  an  isotropic  conductivity  is  used  both  inside  and  outside  the  spot.  In  this 
respect,  our  model  is  like  Parker's  (1974a).  Parker  carried  out  calculations  only  for 
very  shallow  spots,  however,  whereas  the  present  work  has  been  done  for  a  full  range 
of  aspect  ratios. 

In  addition  to  filling  a  gap  in  a  catalog,  the  work  presented  here  has  a  second 
objective:  to  develop  the  simplest  possible  model  of  the  inhibition  concept.  In  the 
absence  of  an  accurate  theory  of  convective  heat  transport,  one  can  make  a  strong 
case  for  simplicity.  The  principal  advantage  of  a  simple  model  is  that  it  allows  a 
thorough  discussion  of  the  dependence  of  the  solution  on  the  basic  spot  parameters. 

The  mathematical  formulation  for  the  present  model  is  given  in  Section  2.  In 
Section  3  some  two-dimensional  models  are  analyzed  in  detail.  The  analysis  leads  to 
some  simple  asymptotic  approximations  which  are  used  to  advantage  in  the  dis¬ 
cussion  of  three-dimensional  models  (Section  4).  The  reader  not  interested  in 
mathematical  details  is  invited  to  skip  to  Section  5,  where  conclusions  from  this  work 
and  other  thermal  models  are  discussed.  In  order  to  make  Section  5  reasonably 
complete,  the  discussions  of  the  solution,  which  would  normally  be  scattered  through 
Sections  2,  3,  and  4,  are  deferred  until  Section  5. 


2.  Formulation 


The  geometry  of  the  model  is  shown  in  Figure  1 .  We  assume  that  the  heat  flux  Q  and 
the  superadiabatic  temperature  gradient  are  linearly  related.  The  thermal  conduc¬ 
tivity  is  assumed  isotropic,  with  a  reduced  value  in  the  volume  Vx  where  the 
inhibition  is  effective.  Thus  we  take 


Q  = 


- K(VT-rk )  in  V2 
-K(l-e)(VT-rk)  in 


Vi, 


(1) 


z 


Fig.  1 .  Model  geometry.  The  volume  Vt  is  the  region  of  low  thermal  conductivity,  and  V2  is  the  region  of 
high  conductivity,  with  S  the  interface  between  them.  The  gas  radiates  into  the  half-space  z  <  0  from  the 
surface  z  =  0,  which  is  divided  into  S,  (bounding  V,)  and  S2  (bounding  V2).  The  heat  flux  is  asymptotically 

equal  to  F ,r>  as  z -*  +oo. 
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where  T  is  the  temperature,  F  is  the  adiabatic  gradient  (assumed  constant),  k  is  a  unit 
vector  in  the  z  -direction,  K  is  the  thermal  conductivity  outside  the  spot,  and  e  is  the 
fractional  reduction  in  conductivity  within  the  spot.  Both  K  and  e  are  taken  to  be 
spatially  constant.  The  steady-state  heat  equation  is  V  •  Q  =  0,  and  this  becomes 

V2r  =  0  inV,andV2.  (2) 

The  temperature  and  the  normal  component  of  heat  flux  must  be  continuous  at  the 
interface  5  and  Vt  and  V2.  Thus 

[71?  =  0 

and 

n  *  (VT-rk)|2  =  (l-e)n  •  (VT  — J\)|i  on  S.  (3) 

At  the  upper  surface,  the  gas  radiates  according  to  the  Stefan -Boltzmann  law,  so  we 
have 

on  S, 
and 

^(~r)  =  o-r4  on  s»,  (4) 


where  cr  is  the  Stefan-Boltzmann  constant.  The  final  condition  is  the  asymptotic  one 
far  beneath  the  spot,  where  the  heat  flux  must  approach  the  solar  flux  F©.  Thus 


dT 

dz 


(5) 


(Here  our  approach  differs  from  that  of  Eschrich  and  Krause  (1977),  whose  imposed 
exact  flux  uniformity  at  the  base  of  their  cylindrical  volume  VV) 

In  the  absence  of  a  spot  (e  -  0)  the  temperature  gradient  is  constant,  and  the  flux 
and  surface  temperature  TSo  are  related  by 

F@  =  <rTtS0.  (6) 


The  temperature  distribution  is  then 
To(z)  =  Tso(l+z/h) , 

where 


h  = 


^i^V‘  -V  -  i  gZjsV* 

vrdz/.-o  vrso  k  ) 


(7) 

(8)  1 

j 


is  the  temperature  scale  height  at  the  surface.  Another  important  scale,  the 
superadiabatic  temperature  scale  height  at  the  surface,  is  defined  by 


trTso 


(9) 
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Returning  to  the  general  case,  we  introduce  Ts o  as  a  temperature  scale  and  L  as  a 
length  scale  (to  be  chosen  later,  but  always  of  the  order  of  the  spot  diameter).  Also, 
we  separate  out  the  behavior  of  T  at  z  *  oo.  Thus  we  introduce 

z  =  L~xz,  V  =  LV, 

and  a  new  dependent  variable  <P ,  defined  by 

r  =  rso(i+|f+j0),  do) 

where 


A=(4  L/H)  (11) 

is  of  the  order  of  twice  the  spot  diameter  divided  by  the  superadiabatic  temperature 
scale  height.  With  this  scaling,  the  problem  becomes: 

VJ0  =  O  in  Vi  and  V2 ,  (12) 

0-*O  as  z->oo,  (13) 

[<P]2i  =  0  on  S,  (14) 

and 


n  •  -  (1  -e)n  •  V0\t -n  •  k  on  S. 

The  upper  surface  boundary  condition  is  now  given  by 

f(l  +Ae#/4)4-  1  +  e 


30 

dz  " 


e(l-e) 


(l+A*#/4)4-l 


on 


on  52 . 


(15) 


(16) 


Because  of  the  boundary  condition  (16),  the  problem  is  nonlinear.  The  solution 
would  require  elaborate  numerical  work,  which  can  hardly  be  justified  in  view  of  the 
other  limitations  of  the  model.  Thus  we  linearize  the  problem,  based  on  the 
assumption  of  small  e.  This  restricts  our  ability  to  make  quantitative  statements 
about  temperature  amplitudes,  but  we  can  still  have  some  confidence  in  the  spatial 
scales  emerging  from  the  calculation.  In  Section  5,  we  discuss  a  crude  but  not 
unreasonable  way  to  restore  some  of  the  nonlinearity  omitted  from  the  calculations. 
(Some  comments  on  other  models  are  in  order  here.  Parker  (1974a)  also  linearized 
the  radiation  condition.  In  his  case,  the  small  parameter  was  the  ratio  of  the  spot 
depth  to  the  temperature  scale  height.  In  the  present  work,  we  have  preferred  to  take 
the  conductivity  contrast  e  to  be  small,  since  we  wish  to  study  deep  spots.  Eschrich 
and  Krause  (1977)  also  linearized  the  radiation  condition.  From  the  point  of  view  of 
mathematical  consistency,  their  linearization  cannot  be  justified,  since  neither  the 
spot  depth  nor  the  conductivity  contrast  is  small  in  their  calculations.  However,  as 
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they  point  out,  the  error  associated  with  their  linearization  is  probably  less  than 
25%.) 

The  linearization  is  carried  out  by  expanding  4>  in  powers  of  e  and  keeping  only 
the  first  term  in  the  expansion.  The  result  is: 


=  0  in  Vj  and  V2 , 

(17) 

<P-*0  as  z-*a o. 

(18) 

[<P]]  =  0,  n  •[$<£]?  = -n  ♦  k  on  S, 

(19) 

d4>  f  1  on  Si 

— -A0-|o  on  52 

(20) 

In  this  formulation,  it  becomes  clear  that  the  character  of  the  solution  depends  only 
on  the  shape  of  the  spot  and  the  parameter  A.  As  we  shall  see  later,  the  case  of  interest 
is  A  » 1. 

Although  the  problem  defined  by  (17)— (20)  is  linear,  it  is  made  somewhat  difficult 
by  the  jump  conditions  on  the  internal  interface  S.  The  basic  technique  we  use  is  first 
to  find  a  harmonic  function  x  which  satisfies  the  jump  condition  (19)  and  the 
condition  (18)  at  oo,  but  not,  in  general,  the  surface  boundary  condition  (20).  Thus  we 
let  x  be  a  function  which  satisfies 

$2*=0  inVjandV2,  x  0 ,  (21) 

[*]?  =  0,  «•[$*]?= -a- It  on  S.  (22) 

To  make  x  unique,  we  also  specify  that 

*=0  on  $!  +  S2.  (23) 


Once  we  have  found  such  a  *,  we  define  a  new  dependent  variable  ^  by 
<P  =  X  +  V 

Then  ^  is  a  solution  of  the  following  boundary  value  problem: 

=  0  in  Vx+V2 , 

0  as  z  oo  , 

and 

a W 

— -A  ¥  =  f  on  f  =0 , 

where 


/= 


dz 

dz 


1 

on 


on  Si 

S2. 


(24) 

(25) 

(26) 

(27) 

(28) 
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The  problem  (25M28)  involves  no  internal  interfaces  and  is  an  ordinary  half -space 
problem  which  can  be  solved  by  the  Fourier  transform.  Thus  we  can  deal  with  any 
spot  geometry  for  which  a  ‘jump’  function  x  can  be  found.  As  we  shall  see,  this  is 
easily  done  for  elliptical  spots  in  the  two-dimensional  case,  and  both  prolate  and 
oblate  spheroidal  spots  in  the  three-dimensional  case. 


3.  Two-Dimensional  Models 

3.1.  Solution  of  half-space  problem 

In  this  section,  we  solve  the  problem  (25)— (28)  for  V'U ,  z),  where  x  is  the 
horizontal  coordinate.  We  choose  the  length  scale  L  so  that  the  spot  volume  Vx 
intersects  the  surface  z  -  0  along  -  1  s  <  1 .  The  dimensional  spot  width  is  then  2L. 
The  problem  is  easily  solved  by  a  Fourier  transform.  With  the  aid  of  the  convolution 
theorem,  one  can  put  the  solution  in  the  form 


^(x 


,z)=  J  G(x-x\z)f( x')dr\ 


(29) 


where 


G(x 


1 ' 


~(ikx  +-|fc|z) 


|fc|  +  A 


dk. 


(30) 


It  is  easy  to  show  (with  the  formulas  in  Chapter  5  of  Abramowitz  and  Stegun,  1964) 
that  the  kernel  G(x,  z)  can  be  expressed  in  terms  of  an  exponential  integral  Ex  : 


where 


and 


GU,7)  =  -Rej^£',(^)|, 


£  =  A  (z+ix). 


(31) 


(32) 


£.(£)=  J  yd/. 


(33) 


An  important  special  case  is  G(jc,  0),  which  is  needed  in  the  calculation  of  the  surface 
temperature.  By  starting  from  (30)  with  z  -  0,  one  can  show  that 


G(x,  0)  =  --g(A|x|), 

1 T 


(34) 


where  g  is  an  auxiliary  function  of  the  exponential  integrals,  defined  by  (Abramowitz 
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and  Stegun,  1964;  p.  232) 

g(M)  =  -cos(w)C/(w)-sin(w)5/(w) ,  (35) 

where 

u 

r  cos  t  —  i 

Ci(w)  =  6 +  ln(w)  +  J  — - — dr  (36) 

o 

and 


o 


(37) 


where  6=0.577  215  664  9  is  Euler’s  constant.  For  small  u,  the  most  useful 
representation  of  g  is  obtained  from  the  expansions 


OO 


C/(M)  =  a  +  ln(«)+  I 

n  *  1 


(-i)V" 

2n(2n)\ 


(38) 


and 


j/(m)=--+  £ 
l  ««0 


(2n  +  1  )(2«  + 1)! ' 


(39) 


For  1  <  u  <  oo,  the  most  useful  representation  is  the  rational  approximation  (Has¬ 
tings,  1955) 


,  v  1  rw8  +  fl1w6  +  a2M4  +  a3K2  +  a4-|  ,  ,  , 

*U)’=?L,  +  i1«6+M4  +  M1  +  *4J  *(“)> 


(40) 


where 


flU  =42.242  855, 


b\  =48.196  927, 


a2  =  302.757  865, 
a3  =  352.018  498, 
Qa  =  21.821  899, 


b7  =  482.485  984, 
63=  1114.978  885, 
64  =  449.690  326, 


(41) 


with  |e(w)|  <  3  x  10~7  for  1  <  u  <  00. 

The  parameter  A  =  4 L/H  is  important  in  the  solution,  so  we  digress  briefly  to 
discuss  typical  values.  From  the  model  of  the  convection  zone  given  by  Spruit 
(1977b),  we  estimate  that  //  —  500  km  at  t  =  1.  The  diameter  of  the  region  of 
inhibition  of  convection  in  our  model  is  2 L,  and  it  is  logical  to  identify  this  with  a 
typical  umbral  diameter  -  say  2L  —  20  000km.  This  gives  A— 80.  There  are 
uncertainties.  In  particular,  H  may  well  be  different  in  other  convection  zone 
models,  since  it  is  determined  by  the  superadiabatic  gradient  in  the  uppermost  layers 
of  the  convection  zone  -  a  particularly  difficult  quantity  to  predict.  Observed  umbral 


24 


.7 


THERMAL  MODELS  OF  SUNSPOTS 


313 


diameters  vary  widely,  so  L  is  also  uncertain.  In  spite  of  the  uncertainties,  it  is  clear 
that  A  is  large,  and  that  is  all  that  is  necessary  for  the  calculations  here.  When  a 
numerical  value  is  required,  we  have  taken  A  =  100  for  purposes  of  illustration. 

It  is  worth  noting  that 


—  [  g(|w|)dH  =  l, 
7T  J 


a  result  that  is  most  easily  proved  from  the  Fourier  transform  of  G(x ,  0).  Thus  for 
large  A,  we  expect 

lim  — g(A|jt-x'|)  =  S(x-x') .  (43) 

A  —ao  TT 

The  formulas  given  here  reduce  the  problem  for  ^  to  quadrature.  To  go  further, 
we  must  choose  a  geometry  for  the  spot  and  find  the  corresponding  jump  function  *. 

3.2.  Elliptical  spots 

We  seek  a  function  *(x,  z)  satisfying  (21)-(23),  when  the  interface  S  is  given  by 

jc2  +  (z/£>)2=  1  ,  (44) 

where  D  is  the  aspect  ratio  of  the  spot.  The  determination  of  a  suitable  \  is 
straightforward,  although  tedious,  and  can  be  carried  out  with  the  aid  of  formulas  for 
the  solution  of  Laplace’s  equation  in  elliptic  coordinates  (Morse  and  Feshbach,  1 953, 
p.  1195-1200).  We  skip  the  lengthy  derivation,  and  only  quote  the  results.  We  begin 
with  shallow  spots  (D  <  1).  We  have 


where 


X(x,z)  = 


m  =  cosh 


z/(Z>  + 1)  in  Vl 
[D/(l-D2)'/2]e~m  sin  6  in  V2 , 


2{1  -D2}' 


with  m  >  0,  and  where 


,  [  [(X  +  {1  -  P2}‘/2)2  +  z2]1/2  -  [U  -  {1  -  D;}l/2)2  +  Z2]"2' 
\  2{!-£>2}1/2 


e  -  cos' 


with  0 ^  0  s  7T.  The  inversion  of  (46)  and  (47)  is  also  useful: 
x  =  (1  ~D2)W2  cosh  m  cos6y 
z  -  ( 1  —  D2)u 2  sinh  m  sin  6 . 


nESSE^'zm 
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3.3.  Surface  temperature 

We  will  show  that  for  large  K,  a  very  simple  approximation  gives  useful  information. 
This  approximation  proves  invaluable  in  the  analysis  of  three-dimensional  spots  in 
Section  4.  In  addition,  our  calculations  will  elucidate  the  nature  of  the  umbral 
boundary,  whose  sharpness  has  been  the  object  of  much  discussion  in  the  literature. 
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Wc  begin  with  the  exact  expression  for  tf'Xx,  0),  obtained  from  (29)  and  (34): 


1P(*,0)  =  --  f  g(A|jr— *'!)/(*') dx\ 

TT  J 

-oc 


(56) 


The  fact  that  A  is  large  suggests  that  we  use  the  delta-function  approximation  (43)  for 
g.  This  yields 

A^(x,0)^-/(x).  (57) 

This  result  is  also  strongly  suggested  by  the  boundary  condition  (27).  The  approxi¬ 
mation  (57)  is  very  attractive,  since  /(x)  is  known  in  closed  form.  However,  it  cannot 
be  uniformly  valid  in  x,  since  f(x)  is  discontinuous  at  x  -±1,  whereas  ip(x,  0)  is 
continuous.  This  is  the  classic  signal  of  a  singular  perturbation  problem,  and  we  must 
use  stretched  coordinates  in  the  vicinity  of  x  =  ±1  in  order  to  remedy  the  defects  in 
(57)  (see,  e.g..  Van  Dyke,  1975,  for  a  general  discussion).  Since  tf'fx,  0)  is  even  in  x, 
we  need  consider  only  x  =  1 .  There,  as  one  can  show,  the  proper  stretched  coordinate 
is 

tj  =*  A(x  —  1) .  (58) 

We  introduce  this  and  V  =  A  (x'  -  1 )  in  (56).  We  then  take  the  limit  A  -»  oo  at  fixed  17. 
After  some  elementary  manipulations  in  which  we  make  use  of  (42),  we  can  put  the 
result  in  the  form 

ini 

lim  A^(x, 0)  =  0(t})  =  — i(/l  +  A)- (— — ^-)sgn(Tj)[  g(u)du,  (59) 
A  -*<30  \  IT  /  J 

v  fixed  0 

where  /*  =  /(l*)  are  the  right  and  left  limits  of  /  at  x  =  1.  The  inner  solution  fi(v) 
blends  smoothly  with  the  outer  solution  A  '/'(x,  0),  since 


lim  /2(rr)  =  -/±.  (60) 

TJ  -+±00 


a  result  easily  established  from  (42)  and  (59).  The  smooth  variation  of  A  V  across  the 
umbral  boundary  is  given  by  (59).  We  see  that  A  'V  varies  from  f  to  f,  for  4»j  ~  1. 
Thus  from  (58)  we  conclude  that  the  thickness  of  this  transition  region  is  of  the  order 
of  A  1 ,  where  A  is  given  by  (1 1).  Thus  the  horizontal  width  of  the  region  of  rapidly 
changing  temperature  near  the  umbral  boundary  is  of  the  order  of  the  superadiabatic 
temperature  scale  height. 

For  computational  purposes,  it  is  highly  desirable  to  blend  (57)  and  (59)  into  a 
uniformly  valid  approximation  (again  see  Van  Dyke,  1975,  for  a  general  discussion). 
This  is  easily  done  and  the  result  is 


A^.=s  I  -f(x)  +  D(ri)+f-,  v<0 

\-f(x)  +  f)(ri)  +  ft,  17  > 0  . 


(61) 


27 


316 


ALFRED  CLARK.  JR 


We  now  have  three  representations  for  <P(x,  0):  the  simplest  approximation  (57), 
the  more  complete  approximation  (61),  and  the  exact  representation  (56).  From 
their  derivation,  we  know  that  (57)  and  (61)  are  asymptotically  valid  for  large  A. 
However,  only  by  computation  and  comparison  with  the  exact  solution  (56)  can  we 
determine  whether  the  approximations  are  numerically  useful  for  A  values  of  interest 
(A  —  100).  For  this  reason  we  have  carried  out  a  detailed  numerical  comparison  of 
(56),  (57),  and  (61).  The  numerical  integration  of  (56)  involves  four  difficulties:  (i)  the 
infinite  range  of  integration;  (ii)  the  large  parameter  A,  requiring  a  fine  grid;  (iii)  the 
logarithmic  singularity  of  g  at  x  =  x';  (iv)  the  discontinuities  in  f  at  x  =  ±  1 .  A  detailed 
description  of  the  numerical  scheme  actually  used  would  take  many  pages.  Since 
there  are  no  major  new  features  in  the  scheme,  we  forgo  the  description. 

The  relation  between  the  surface  temperature  7s(x)  and  tf'tx,  0)  is  easily  obtained 
from  Equations  (10),  (23),  and  (24): 


Ts(x)-Ts o  v  A^(x,0) 

- ~ - -  dTs(x)  = - - - 

eTS0  4 


(62) 


Here  ATs(x)  is  the  fractional  change  in  surface  temperature  per  unit  fractional 


TABLE  l 

The  change  in  surface  temperature  ATs(x)  for  A  =  100  and  D  =  5 


Values  of  arsU) 

Exact, 

Approximate, 

Approximate, 

X 

based  on  (56) 

based  on  (61) 

based  on  (57) 

0.0 

-0.2079 

-0.2075 

-0.2083 

0.2 

-0.2078 

-0.2073 

-0.2083 

0.4 

-0.2076 

-0.2070 

-0.2083 

0.6 

-0.2070 

-0.2063 

-0.2083 

0.8 

-0.2051 

-0.2044 

-0.2083 

1.0 

-0.0836 

-0.0833 

undefined 

1.2 

0.0362 

0.0357 

0.0397 

1.4 

0.0362 

0.0358 

0.0378 

1.6 

0.0350 

0.0346 

0.0359 

1.8 

0.0335 

0.0331 

0.0341 

2.0 

0.0320 

00316 

0.0324 

2.2 

0.0304 

0.0301 

0.0307 

2.4 

0.0289 

0.0286 

0.0292 

2.6 

0.0275 

0.0272 

0.0277 

2.8 

0.0261 

0.0258 

0.0262 

3.0 

0.0248 

0.0245 

0.0249 

3.2 

0.0235 

0.0232 

0.0236 

3.4 

0.0223 

0.0221 

0.0224 

3.6 

0.0212 

0.0209 

0.0212 

3.8 

0.0201 

0.0199 

0.0202 

4.0 

0.0191 

0.0189 

0.0191 

28 
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TABLE  II 

The  change  in  surface  temperature  drs(jc)  for  x  near  1,  and  for  A  *  100  and  D  =  5 


Values  of  4TsU) 

Exact, 

Approximate, 

Approximate, 

X 

based  on  (56) 

based  on  (61) 

based  on  (57) 

0.70 

-0.2064 

-0.2057 

-0.2083 

0.80 

-0.2051 

-0.2044 

-0.2083 

0.85 

-0.2038 

-0.2031 

-0.2083 

0.90 

-0.2013 

-0.2005 

-0.2083 

0.91 

-0.2004 

-0.1997 

-0.2083 

0.92 

-0.1994 

-0.1987 

-0.2083 

0.93 

-0.1981 

-0.1974 

-0.2083 

0.94 

-0.1964 

-0.1957 

-0.2083 

0.95 

-0.1941 

-0.1934 

-0.2083 

0.96 

-0.1909 

-0.1901 

-0.2083 

0.97 

-0.1859 

-0.1851 

-0.2083 

0.98 

-0.1774 

-0.1766 

-0.2083 

0.99 

-0.1597 

-0.1589 

-0.2083 

1.00 

-0.0836 

-0.0833 

undefined 

1.01 

-0.0075 

-0.0079 

0.0416 

1.02 

0.0101 

0.0097 

0.0415 

1.03 

0.0185 

0.0181 

0.0414 

1.04 

0.0234 

0.0230 

0.0413 

1.05 

0.0266 

0.0262 

0.0412 

1.06 

0.0288 

0.0284 

0.0411 

1.07 

0.0304 

0.0300 

0.0410 

1.08 

0.0316 

0.0312 

0.0409 

1.09 

0.0326 

0.0321 

0.0408 

1.10 

0.0333 

0.0329 

0.0407 

1.15 

0.0354 

0.0349 

0.0402 

1.20 

0.0362 

0.0357 

0.0397 

1.30 

0.0365 

0.0361 

0.0387 

change  in  thermal  conductivity.  Tables  1  and  II  give  a  numerical  comparison  of  the 
two  approximations  and  the  exact  solution  for  the  parameter  values  A  =  100  and 
D  =  5.  It  can  be  seen  that  the  uniformly  valid  approximation  (61)  is  in  excellent 
agreement  with  the  exact  solution.  Even  the  simple  approximation  (57)  is  very  good 
except  in  the  immediate  vicinity  of  x  *  1 .  Calculations  for  other  parameter  values 
(A  =  50  and  100;  D  =  0.5, 1.0,  2.0,  5.0,  and  10.0)  give  similar  results.  Outside  of  the 
thermal  transition  region  of  the  umbral  boundary,  the  simple  approximation  (57)  is 
quite  adequate.  We  will  take  advantage  of  this  in  our  discussion  of  three-dimensional 
spots  in  Section  4. 

Because  the  sharp  umbral  edge  is  of  some  interest,  the  solution  near  x  =  1  (both 
exact  and  approximate)  is  shown  in  Figure  2,  again  for  the  parameter  values  A  =  100, 
D  =  5.  As  predicted  by  the  asymptotic  theory,  the  width  of  the  region  of  sharp 
temperature  change  is  of  the  order  of  A  ~ 1 . 
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Fig.  2.  Surface  temperature  distribution  in  the  vicinity  of  the  umbral  edge  (x  =  l).  The  center  of  the 
umbra  is  at  x  =  0.  The  aspect  ratio  of  the  spot  is  D  =  5.  The  value  of  A  is  100,  corresponding  to  a  spot 
diameter  of  50  times  the  superadiabatic  temperature  scale  height.  The  solid  line  is  the  exact  solution  (56) 
and  the  uniformly  valid  approximation  (61),  which  are  indistinguishable  on  this  scale.  The  dashed  line  is 

the  simple  approximation  (57). 


3,4.  Variation  of  temperature  with  depth 


Although  the  surface  temperature  distribution  is  the  principal  observable  quantity, 
the  distribution  of  temperature  below  the  visible  surface  is  also  of  interest.  In  this 
section,  we  obtain  some  important  qualitative  results  about  the  depth  dependence  of 
the  temperature. 

The  starting  point  is  Equation  (10)  for  the  temperature,  Equation  (24)  for  <Py 
Equations  (29)— (33)  for  and  Equations  (45)—(55)  for  *U,  z)  and  /(*).  The 
unperturbed  temperature  r0(z )  is  given  by  (7).  A  useful  measure  of  the  temperature 
distribution  is  AT( x,  z),  defined  to  be  the  fractional  change  in  temperature  per  unit 
fractional  change  in  conductivity: 


AT(x  ,z)« 


ru,  z)-r0(z)  a(x  +  ^) 
eTo(z)  4(1  -bAz) 


(63) 


For  z  =0,  this  reduces  to  ATs(x)  -  A^(jc,  0)/4,  discussed  in  the  previous  section. 
Away  from  the  surface,  however,  the  term  x  dominates.  To  show  this,  we  use  the 
expression  (31)  for  G,  and  observe  that,  for  A  »  1,  and  Az  »  1,  the  complex  argument 
f  -  A[z  +  /(jt  -*')]  is  always  large.  Then  we  can  use  an  asymptotic  expansion  for 
Ei(f)  (Abramowitz  and  Stegun,  1964;  p.  231)  to  write 
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Substitution  of  this  into  Equation  (29)  for  IP  yields 

<30 

V(x,z)  =  ~  [  -  fiX  ),^  i  +  0(\-2)  for  A*  » 1 .  (65) 

irk  J  (x  -  x  )  +  z 

-00 

Thus  V  is  small  of  order  A-1  provided  Az  »  1.  Since  x  is  of  order  one  with  respect  to 
A,  we  conclude  that 

for  Az»l.  (66) 

4  z 

Since  x  >  0  for  z  >  0,  we  get  the  important  conclusion  that  the  temperature  pertur¬ 
bation  is  positive  away  from  the  upper  boundary.  Since  ^  is  0(A  " 1 )  even  on  z  -  0,  we 
have  V  =  0(k  ~*)  everywhere,  and  x  will  dominate  except  where  x  =  0.  It  is  easy  to 
show  that  x  ^  0(k  -I)  only  for  Az  1 .  Thus  the  cooling  of  a  spot  is  purely  a  surface 
phenomenon ;  only  within  one  or  two  superadiabatic  temperature  scale  heights  of  the 
visible  surface  is  the  spot  cooler  than  the  normal  atmosphere.  The  remainder  of  the 
interior  of  the  spot  is  hotter  than  the  normal  atmosphere.  A  more  quantitative 
discussion  of  this  effect  is  given  in  Section  5.2  for  three-dimensional  spots. 


4.  Three-Dimensional  Spots 


We  choose  half  spheroids  for  the  geometry  since  they  are  sufficiently  simple  to  allow 
an  analytical  solution,  but  sufficiently  general  to  include  deep  and  shallow  spots.  We 
choose  our  basic  length  scale  L  to  be  the  radius  of  the  horizontal  circular  section  of 
the  spheroid.  The  equation  of  the  interface  5  is  then 


x2+y2+j?= 1  • 


(67) 


The  discussion  in  Section  3  shows  that  all  essential  information  can  be  obtained 
from  the  jump  function  x  by  means  of  simple  approximations.  Thus  we  limit  our 
calculations  to  the  determination  of  x «  As  in  the  two-dimensional  case,  the  cal¬ 
culation  is  lengthy  but  straightforward.  The  necessary  formulas  are  given  by  Morse 
and  Feshbach  (1953,  pp.  1284-1294).  As  before,  we  skip  the  derivation  and  only 
quote  the  results.  We  begin  with  deep  spots  (prolate  spheroids;  D  >  1).  Then  x  is 
most  conveniently  expressed  in  terms  of 


((r,z)  = 


{[2  +  (D2  -  1),/2]2  +  r2}1/2  +  {[z  -  (D2  - 1  ),/2]2  +  r2},/2 
2[DT-V^ri 


(68) 


where  r  =  (jc2  +  y2)'/2 


is  the  cylindrical  radial  coordinate.  The  interface  S  is  given 


by 


D 

(D2-  1)1/2’ 


(69) 
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The  special  case  of  spherical  spots  (D  =  1)  can  also  be  treated.  One  gets 


X(r,z)  = 


in  Vt 
3 

l(TTP)575  in  V2’ 


(78) 


and  the  boundary  function  /  is 
2 

3  > 


/(')  = 


*  r<\ 


-p, 


(79) 


Consider  now  the  surface  temperature.  We  use  the  approximation  (57)  and  the 
definition  (62)  to  get 


ATs(r)  = 


Ts(r)-Tso  AV>(r,0)  f(r) 

eTS0  4  4 


(80) 


In  this  approximation,  we  do  not  resolve  the  sharp  but  smooth  transition  at  the 
umbral  boundary,  but  we  do  get  good  estimates  of  (i)  the  umbral  temperature,  (ii)  the 
maximum  temperature  in  the  bright  ring,  and  (iii)  the  horizontal  scale  for  the  fall-off 
of  the  temperature  in  the  bright  ring  as  we  move  away  from  the  spot.  In  addition,  this 
approximation  satisfies  the  global  conservation  of  energy.  Let  AF  = 
4<rTlo(Ts-TSQ)  be  the  (linearized)  perturbation  in  the  radiative  flux  at  the  free 
surface.  Then  the  fractional  flux  perturbation  is 

AF  4 (Ts-TSo)  v  ,r#/  £  /otv 

— r-  - - - lP(r,  0)  ~  -ef.  (81) 

<7*  SO  *S0 

Conservation  of  energy  requires  that  the  integral  over  Si  +  52  of  AF  vanish.  It  is 
straightforward  to  show  this,  by  using  the  definition  (28)  for  /,  the  jump  conditions 
(22)  for  x*  the  Equation  (21)  for  and  the  divergence  theorem. 

Finally,  it  is  worth  noting  that  within  the  framework  of  the  large  A  approximation 
used  here,  the  character  of  the  solution  depends  only  on  the  aspect  ratio  D. 


5.  Discussion 

The  discussions  based  on  the  calculations  of  Sections  2,  3,  and  4  all  have  been 
deferred  to  this  section,  in  order  to  present  them  in  a  coherent  and  unified  way#. 
Specific  points  of  interest  are  discussed  in  Sections  5.1  (the  sharpness  of  the  umbral 
boundary),  5.2  (the  internal  temperature  of  spots),  5.3  (bright  rings  and  penumbrae), 
and  5.4  (spot  depth).  Concluding  remarks  and  a  summary  are  presented  in  5.5. 

5.1.  The  sharpness  of  the  umbral  boundary 

Several  authors  (e.g.,  Cowling,  1953;  Parker,  1974a)  have  suggested  that  the 
observed  sharp  umbral  boundary  is  inconsistent  with  models  based  on  the  diffusion 
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equation,  unless  the  depth  of  the  region  of  thermal  inhibition  is  very  small.  However, 
very  shallow  spots,  in  the  inhibition  theory,  produce  intense  bright  rings  (Parker, 
1974a)  and  are  therefore  inconsistent  with  observation  (see  Section  5.4  below  for 
further  discussion  of  spot  depth).  Cowling  (1976a,  b)  has  suggested  that  it  is  not 
necessary  for  the  depth  to  be  small,  provided  that  the  horizontal  thermal  conductivity 
in  the  region  of  strong  field  is  greatly  reduced.  Even  this  constraint  is  not  necessary, 
as  the  present  model  shows.  Of  course  the  horizontal  conductivity  in  the  spot  may  be 
small.  The  point  is  that  the  sharp  umbral  boundary  does  not  require  it.  In  fact,  it  is 
difficult  not  to  get  a  sharp  umbral  boundary  out  of  thermal  models  -  the  models  of 
Eschrich  and  Krause  (1977),  Spruit  (1977a)  and  the  present  work  all  show  sharp 
umbral  boundaries.  We  agree  completely  with  Spruit  (1977a)  who  says:  “To 
maintain  a  sharp  transition  in  temperature,  a  sharply  bounded  and  strong  inhomo¬ 
geneity  in  the  diffusion  coefficient  is  sufficient.**  We  would  add  to  this  that  the 
conductivity  inhomogeneity  must  be  present  at  the  visible  surface,  since  it  is  really 
the  combination  of  the  inhomogeneity,  the  radiative  boundary  condition,  and  the 
small  superadiabatic  temperature  scale  height  (compared  to  the  spot  size)  which 
gives  the  sharp  umbral  boundary. 

Neither  Spruit  nor  Eschrich  and  Krause  identify  the  thickness  of  the  umbral 
temperature  transition  in  terms  of  other  scales.  As  our  analysis  in  Section  3.3  shows, 
this  horizontal  transition  region  has  a  thickness  equal  to  the  superadiabatic  tempera¬ 
ture  scale  height  at  the  surface  (defined  by  (9)),  provided  that  the  spot  dimensions 
greatly  exceed  this  scale. 

It  should  be  emphasized  that,  even  though  thermal  models  predict  the  sharp 
umbral  boundary,  they  do  not  explain  it.  What  they  show  is  that  if  there  is  a  sharp 
change  in  the  thermal  conductivity,  then  there  will  be  a  sharp  change  in  temperature 
at  the  visible  surface.  In  the  inhibition  theory,  a  sharp  change  in  thermal  conductivity 
is  caused  by  a  sharp  change  in  magnetic  field  strength.  Interestingly  enough,  the 
theory  of  cooling  by  Alfven  waves,  as  developed  by  Parker  (1974b),  also  requires  a 
sharp  change  in  magnetic  field  strength  to  get  a  sharp  change  in  temperature.  Thus 
while  both  theories  can  accommodate  a  sharp  umbral  boundary,  there  remains 
the  fundamental  question  of  why  the  magnetic  flux  tube  of  the  spot  has  a  sharp 
boundary. 

5.2.  The  internal  temperature  of  sunspots 

As  Parker  (1974a)  has  shown  for  shallow  spots,  and  as  the  present  work  (Section  3.4, 
Section  4)  shows  for  both  deep  and  shallow  spots,  the  subsurface  interior  of  a  spot  is 
hotter  than  the  normal  atmosphere  at  the  same  level.  In  fact,  in  a  thermal  model  of  a 
spot,  the  cool  region  extends  downward  only  a  distance  of  the  order  of  the 
superadiabatic  temperature  scale  height.  It  is  of  interest  that  considerations  of 
mechanical  equilibrium  also  suggest  that  sunspot  cooling  is  cohfined  to  the  surface 
layers  (Weiss,  1964,  1969). 

Consider  now  an  estimate  of  the  magnitude  of  this  subsurface  heating.  We  begin 
with  the  formula  (10),  from  which  we  compute  the  dimensional  vertical  temperature 
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gradient  in  the  spot  volume  (Vi): 


(82) 


At  distances  of  several  superadiabatic  temperature  scale  heights  or  more  below  the 
visible  surface,  0  —  *  (as  discussed  in  Section  3.4).  For  the  spheroidal  spots  of 
Section  4,  the  gradient  in  Vx  is  constant,  and  we  have 


i 

3 


for 

D<  1 

for 

D  =  1 

for 

D>  1  , 

(83) 


where 


(84) 


Here  the  spot  volume  Vx  is  one -half  of  a  spheroid,  and  D  is  the  ratio  of  the  depth  to 
the  radius  of  the  spheroid.  To  complete  the  estimate,  we  need  a  value  of  e.  We  can 
relate  e  to  the  umbral  temperature  Tu.  From  (10)  evaluated  at  2=0,  with  the 
approximation  (57)  for  and  the  definition  (28)  for  /,  we  get 


We  may  solve  (85)  for  e.  Then  by  substitution  into  (82)  we  get 


d71  Ts  o 
T7L  =-n-(l+“). 


dz  I 


H 


(86) 


where 

a=4(l-ru/rso)(^^).  (87) 

Then  a  is  the  fractional  excess  temperature  gradient  in  the  spot.  For  illustrative 
purposes,  we  choose  Tv  =4000K  and  rso  =  5800K,  and  then  calculate  a  as  a 
function  of  the  spot  aspect  ratio  D.  Some  typical  values  of  a(D)  are  a  (0.5)  =  1.384, 
a(1.0)  =  0.621,  a(2.0)  =  0.261,  and  a(5.0)  =  0.073.  The  numbers  cannot  be  taken 
too  seriously  since  the  value  of  e  implied  by  (85)  is  not  small,  in  violation  of  the 
mathematical  basis  of  the  calculation.  Nevertheless  the  numbers  suggest  that  the 
effect  may  be  particularly  significant  for  shallow  spots.  Since  neither  Spruit  (1977a) 
nor  Eschrich  and  Krause  (1977)  discuss  the  temperature  variation  beneath  the 
surface,  it  is  not  known  to  what  extent  the  present  results  are  model-dependent. 
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Elementary  physical  reasoning  about  heat  flow  around  obstacles  suggests  strongly 
that  some  subsurface  temperature  excess  in  the  spot  region  is  inevitable.  With  a 
conductivity  which  increases  rapidly  with  depth,  such  as  that  used  by  Spruit  (1977a) 
the  excess  temperature  may  be  smaller  than  estimated  here. 

In  the  traditional  version  of  the  inhibition  theory,  one  has  usually  assumed  that  the 
coolness  of  the  spot  extends  over  the  spatial  region  in  which  inhibition  is  important. 
A  very  attractive  feature  of  that  assumption  is  that  it  leads  to  a  mechanism  for  field 
concentration,  first  suggested  by  Parker  (1955).  In  simplest  terms,  the  mechanism  is 
the  following:  increasing  the  field  strength  gives  greater  inhibition  of  convection, 
causing  the  spot  to  become  cooler.  As  the  spot  cools,  its  pressure  drops,  and  it  must 
contract  to  preserve  mechanical  equilibrium.  The  contraction  increases  the  field 
strength  further,  and  the  process  continues  until  some  kind  of  equilibrium  is  reached. 
(See  Galloway  et  ai  (1977)  for  a  recent  discussion  of  this  and  other  field-concen¬ 
tration  mechanisms.)  As  Parker  (1974a,  1976)  has  emphasized,  and  as  the  present 
work  shows,  this  mechanism  cannot  work  in  the  inhibition  theory,  because  the  spot 
region  is  not  cooled  (except  very  near  the  visible  surface)  but  is  heated  by  the 
inhibition  process.  Parker  (1974a)  concludes  that  the  heating  would  in  fact  disperse 
the  field,  and  that  this  is  a  major  objection  to  the  inhibition  theory.  In  the  absence  of  a 
quantitative  theory  of  spot  formation,  it  is  difficult  to  draw  such  a  definite  conclusion. 
We  prefer  the  following  more  conservative  conclusion:  If  the  inhibition  theory  with 
increased  subsurface  temperatures  is  correct,  then  mechanical  equilibrium  (internal 
gas  pressure  plus  magnetic  pressure  equals  external  gas  pressure)  requires  that  at  any 
given  level  the  spot  have  a  lower  density  than  the  surrounding  atmosphere.  Thus 
some  material  must  be  expelled  from  a  flux  tube  during  the  concentration  phase  of  its 
history.  Alternatively,  as  suggested  by  Meyer  et  ai  (1974),  the  mechanical  balance 
may  be  partially  dynamic  if  there  are  significant  large  scale  flows  outside  the  flux 
tube.  In  any  case,  the  above  values  of  a  suggest  that  the  mechanical  balance  problem 
is  more  severe  for  shallow  spots  than  deep  ones. 

A  digression  on  cooling  by  Alfven  waves  is  useful  here.  In  that  theory,  it  is 
perfectly  possible  for  the  cooling  to  extend  over  a  sizable  volume,  rather  than  being  a 
surface  effect  as  in  the  inhibition  theory.  Such  a  volume  cooling  would  make  sunspot 
formation  much  easier  to  understand.  Cowling  (1976a;  see  also  Parker,  1977)  has 
argued,  on  general  thermodynamic  grounds,  that  it  is  very  difficult  to  transform  all  of 
the  missing  heat  flux  into  Alfven  waves.  However,  that  may  not  be  necessary.  As  the 
calculations  of  this  section  show,  deep  spots,  in  the  inhibition  theory,  have  only 
slightly  warm  interiors.  It  is  possible  that  a  modest  flux  of  Alfven  wave  energy,  in 
combination  with  the  inhibition  of  convection,  can  produce  a  deep  spot  with  a 
substantial  cool  volume.  Such  hybrid  theories,  as  exemplified  by  the  work  of  Mullan 
(1974),  deserve  more  attention  than  they  have  received. 

5.3.  Bright  rings  and  penumbrae 

A  point  of  primary  interest  in  any  sunspot  theory  is  the  disposition  of  the  'missing 
flux'.  As  many  authors  have  pointed  out,  the  inhibition  theory  requires  a  bright  ring 
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around  a  spot,  in  which  the  diverted  heat  flux  is  radiated  from  the  visible  surface.  The 
amplitude  and  scale  of  this  bright  ring  are  of  great  interest,  since  they  are,  in 
principle,  observable.  In  this  section,  we  will  summarize  the  predictions  of  the 
present  model,  and  compare  them  with  those  of  Parker  (1974a),  Eschrich  and 
Krause  (1977),  and  Spruit  (1977a).  As  we  shall  see,  the  comparison  with  observation 
is  complex  and  inconclusive. 

We  begin  with  a  discussion  of  the  spatial  scale  of  the  bright  ring.  In  calculating  the 
surface  temperature,  we  use  the  large  A  approximation  (57).  Then  we  get  from  (10), 
with  z=0, 

7’s(r)=rso[l-j/(r)].  (88) 


where  f(r)  is  given  by  (71)  for  deep  spots,  by  (75)  for  shallow  spots,  and  by  (79)  for 
hemispherical  spots.  Figure  3  shows  a  plot  of  the  normalized  excess  temperature  in 
the  bright  ring. 


ATR(r)  = 


Ts(r)-  Tso 

Ts(l)-TSo' 


(89) 


as  a  function  of  r,  for  D  =  0.2,  1.0,  and  5.0.  For  large  r,  ATR  falls  off  like  l/r  \  The 
width  of  the  region  in  which  ATR  is  appreciable  increases  with  increasing  spot  depth. 
A  more  precise  statement  can  be  made  by  computing  the  radius  rw  at  which 
ATR(rw)  =  0A.  The  results,  given  in  Table  III,  show  that 


rw-\-D 


(90) 


for  a  wide  range  of  values  of  D.  Thus  the  width  of  the  bright  ring  is  comparable  with  the 
depth  of  the  spot. 


Fig.  3.  Temperature  excess  J  TR  as  a  function  of  radius  r  in  bright  ring,  normalized  to  unity  at  the  umbral 
boundary  (r  =  1 ).  The  three  curves  are  for  a  shallow  spot  {D  -  0.2),  a  hemi-spherical  spot  {D  =  1 .0),  and  a 

deep  spot  { D  ~  5.0). 
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TABLE  III 

The  width  rw  of  the  bright  ring  as  a  function  of 
spot  aspect  ratio  D 


D 

(rw  -D/D 

D 

{rw-\)/D 

0.2 

1.28 

1.5 

1.10 

0.4 

1.26 

2.0 

1.05 

0.6 

1.22 

3.0 

0.99 

0.8 

1.18 

4.0 

0.95 

1.0 

1.15 

5.0 

0.92 

Consider  now  the  amplitude  of  the  bright  ring.  Here,  unfortunately,  the  linear 
theory  breaks  down.  One  cannot  produce  an  umbral  temperature  of  4000  K  with  a 
small  e.  In  order  to  discuss  amplitudes,  we  use  the  length-scale  information  obtained 
above  from  the  linear  theory,  and  we  use  a  simple  model  which  incorporates  the 
correct  nonlinear  radiative  boundary  condition  on  z  =  0.  We  do  this  by  assuming  a 
functional  form  for  Ts(r).  The  peak  temperature  excess  in  the  bright  ring,  ATP~ 
pTSo ,  is  left  as  a  free  parameter,  to  be  determined  by  the  global  energy  balance.  Thus 
we  assume 


Ts(r)  = 


r  c.. 


rCl 


rso’  r  >  1 1  +  1  +  ((r—l )/  tv  )3 1 


(91) 


We  choose  w  so  that  the  temperature  excess  has  dropped  to  0.1  ATP  for  r  =  1  +  D  (in 
accordance  with  (90)).  This  gives 

w  =  D/  9W\  (92) 


We  impose  the  global  energy  balance: 


lim  |  (T7t T\j  4*  (T7T  2  T\r  dr  cnr T^qR  ^  |  =  0  . 


(93) 


The  substitution  of  (91)  into  (93)  yields 

Tjj'4 
Iso' 

+  20w2(4Jx  +  6/3/2  +  4 p2Jy  +  /?3/4) , 

where 


1  -  (  — )  4  =  2/3  w  (4 /,  +  6  pl2  +  4/3  :/3  +  /3  V4 )  + 

'  *  SO* 


/  _  f  djr  r  _  f  xdx 

"  J  (1  +  jrV*  i  (1+*V 


(94) 


(95) 
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By  residue  theory,  one  shows  easily  that 

v3 '3  9  27  243' 

and 


(Jl ,  72,  j},  /4)  =  .  — )  . 

>/3'3  9  27  243' 


(96) 


(97) 


From  (92),  (94),  (96),  and  (97),  we  obtain  a  single  numerical  relation  between 
( Tu/Tso ),  P  =  ATp/Tso,  and  D.  Figure  4  shows  a  plot  of  versus  D,  for  several 
values  of  Tv,  when  Tsu  =  58 00  K.  The  value  of  J  Tp  is  insensitive  to  the  umbral 
temperature  but  very  sensitive  to  the  aspect  ratio  D.  It  is  clear  that  a  shallow  spot 
(D<  1)  in  the  inhibition  theory  has  an  intense  bright  ring. 


Fig.  4.  Peak  excess  Temperature  (above  photospheric)  in  the  bright  ring  as  a  function  of  spot  aspect 
ratio  £>.  for  a  photospheric  temperature  of  5800  K.  and  an  umbral  temperature  of  4500  K  (curve  1), 
4000  K  (curve  2),  and  3500  K  (curve  3). 
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IWe  compare  these  results  with  those  of  other  models.  Parker  (1974a)  found 

intense,  narrow  bright  rings,  the  width  being  comparable  with  the  temperature  scale 
height.  This  is  not  inconsistent  with  the  present  results,  since  Parker  assumed  that  the 
spot  depth  is  much  less  than  the  temperature  scale  height.  The  present  work  requires 
i  the  opposite  inequality.  Although  Eschrich  and  Krause  (1977)  do  not  identify  the 

|  width  of  the  bright  ring,  nor  do  they  give  any  numbers  for  it,  their  results  presented 

graphically  appear  to  be  consistent  with  the  estimate  (90)  for  the  width  of  the  bright 
ring.  Their  results  for  the  amplitude  of  the  bright  ring  are  somewhat  smaller  than  but 
quite  comparable  to  our  results.  Spruit  s  (1977a)  calculations  show  much  weaker 
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bright  rings.  He  finds,  for  example,  a  bright  ring  with  a  temperature  excess  of  about 
20  K  for  a  cylindrical  spot  of  radius  5000  km  and  depth  1000  km.  As  one  can  show, 
this  implies  that  the  flux  missing  from  the  umbra  is  spread  over  an  area  with  a 
diameter  of  the  order  of  500  000  km.  Spruit  attributes  this  rapid  spreading  to  the 
strong  increase  with  depth  of  the  conductivity  in  his  model. 

Consider  now  the  problem  of  comparing  bright  ring  models  with  observation.  In 
principle,  the  amplitude  and  scale  of  the  bright  ring  are  observable.  In  practice,  the 
situation  is  greatly  complicated  by  the  presence  of  a  penumbra  in  a  real  spot.  For 
purposes  of  the  present  discussion,  we  will  adopt  the  standard  view  that  the  umbra  is 
a  measure  of  the  basic,  vertical,  subsurface  flux  tube,  whereas  the  penumbra  is  a 
vertically  thin  region  in  the  layers  above  r  -  1  where  the  spreading  magnetic  field  is 
predominantly  horizontal.  The  principal  observed  fact  is  that  the  penumbra  is  darker 
than  the  normal  photosphere  whereas  the  models  considered  here  require  the 
greatest  temperature  excess  to  be  seen  just  outside  the  umbra. 

There  are  (at  least)  two  qualitative  ways  out  of  this  dilemma.  First,  one  can  suppose 
that  the  excess  flux  coming  up  just  beneath  the  penumbra  is  converted  into 
mechanical  energy,  which  propagates  away.  Both  the  Evershed  effect  and  penumbral 
waves  are  candidates  suggested  by  observation.  This  idea  cannot  be  tested  obser- 
vationally  until  there  are  detailed  predictions  about  the  nature  of  the  mechanical 
flux,  and,  most  important,  about  where  the  mechanical  flux  is  dissipated  and  becomes 
visible  as  excess  radiative  flux. 

A  second  possibility,  considered  quantitatively  by  Spruit  (1977a),  is  that  the 
horizontal  magnetic  field  of  the  penumbra  also  inhibits  convection.  In  a  thermal 
model,  this  would  appear  as  an  additional  thin  region  of  lowered  conductivity.  The 
entire  region  of  lowered  conductivity  would  then  have  the  appearance  of  an  inverted 
top  hat.  The  difficulty  with  this  picture  is  that,  as  Parker  (1974a),  Eschrich  and 
Krause  (1977),  and  we  have  shown,  thin  regions  of  inhibition  produce  intense  bright 
rings.  In  this  case,  the  intense  bright  ring  would  appear  just  outside  the  penumbra. 
Whether  the  rings  would  be  too  bright  to  be  consistent  with  observations  cannot  be 
said  until  detailed  model  calculations  have  been  carried  out.  If  the  effective  thermal 
conductivity  of  the  convection  zone  increases  as  rapidly  with  depth  as  Spruit  (1977a) 
has  suggested,  then  as  his  calculations  show,  the  bright  ring  outside  a  penumbral 
model  of  this  sort  will  indeed  be  weak. 

5.4.  The  depth  of  the  inhibition  region 

As  the  calculations  in  Section  5.3  show,  the  thermal  inhibition  theory,  in  its  present 
simple  form,  gives  a  relation  (shown  in  Figure  4)  between  umbral  temperature,  bright 
ring  peak  temperature,  and  spot  aspect  ratio.  Because  of  the  masking  effect  of  the 
penumbra  discussed  above,  we  cannot  draw  any  definite  conclusions  from  the 
present  model  about  the  depth  of  the  inhibition  region.  However,  two  of  our  results 
suggest  that  a  shallow  spot  makes  greater  demands  on  the  theory  than  a  deep  spot. 
First,  there  is  the  peak  bright  ring  temperature  (Figure  4).  The  excess  flux  near  the 
umbral  boundary  becomes  very  large  as  D  drops  below  1 .  On  the  other  hand  a  value 
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of,  say,  D  =  5  poses  only  a  small  theoretical  problem,  since  the  excess  flux  cor¬ 
responds  to  a  maximum  excess  temperature  of  the  order  of  50  K.  The  second  result 
suggesting  deep  spots  concerns  the  increased  internal  spot  temperature  discussed  in 
Section  5.2.  The  increase  is  very  large  for  shallow  spots,  but  is  quite  small  for  D  =  5. 

We  may  ask  whether  a  value  D  =  5  or  greater  is  reasonable  from  the  point  of  view 
of  dynamics.  Weiss  (1964)  has  tabulated,  as  a  function  of  depth,  the  magnetic  field 
strength  necessary  to  inhibit  convection  (the  calculation  being  based  on  the  idea  that 
the  magnetic  field  will  seriously  interfere  with  the  convection  when  the  magnetic 
energy  density  is  comparable  with  the  kinetic  energy  density  of  the  convection). 
Weiss  finds  that  a  field  of  less  than  6000  G  will  interfere  with  convection  throughout 
the  convection  zone.  Thus  a  modest  contraction  with  depth  of  a  flux  tube  with  a 
surface  field  of  3000  G  is  sufficient  to  give  a  deep  region  of  inhibition. 

5.5.  Concluding  remarks 

We  summarize  here  the  most  important  features  of  the  present  model: 

(1)  The  umbral  boundary  is  sharp  even  for  deep  spots.  The  horizontal  scale  for 
temperature  variations  across  the  umbral  boundary  is  the  superadiabatic  tempera¬ 
ture  scale  height  at  the  surface.  Since  the  models  of  Eschrich  and  Krause  (1977)  and 
Spruit  (1977a)  also  show  a  sharp  umbral  boundary,  it  seems  safe  to  conclude  that  it  is 
a  general  consequence  of  the  inhibition  theory. 

(2)  The  missing  umbral  flux  shows  up  in  the  model  as  a  bright  ring  around  the 
umbra.  Since  this  is  a  consequence  of  flux  conservation,  it  is  a  general  feature  of  the 
inhibition  theory.  The  intensity  of  the  bright  ring,  however,  is  highly  model  depen¬ 
dent.  In  our  model,  in  the  model  of  Eschrich  and  Krause  (1977),  and  in  Parker’s 
model  (1974a),  shallow  spots  give  intense  bright  rings.  Spruit  ( 1 977a),  however,  finds 
weak  bright  rings  even  for  shallow  spots.  He  attributes  this  to  the  strong  horizontal 
spreading  of  flux  disturbances  associated  wbh  the  depth-dependent  conductivity  in 
his  models.  All  of  the  models  agree  in  predicting  weak  bright  rings  for  sufficiently 
deep  spots. 

(3)  Only  a  thin  surface  layer  of  the  spot  is  cool.  At  a  depth  of  several  super¬ 
adiabatic  scale  heights  or  more,  the  material  directly  beneath  the  visible  umbra  is 
hotter  than  the  surrounding  normal  atmosphere  at  the  same  depth.  Although 
Eschrich  and  Krause  ( 1 977)  do  not  discuss  this  point,  Spruit  (private  communication) 
finds  similar  results,  and  it  is  probable  that  the  subsurface  heating  is  a  general  feature 
of  the  inhibition  theory. 

In  spite  of  the  different  descriptions  of  turbulent  transport  used  by  Eschrich  and 
Krause  (1977),  by  Spruit  (1977a),  and  in  the  present  work,  all  of  the  models  yield 
surface  temperature  distributions  resembling  sunspots.  This  provides  support  for  the 
inhibition  theory  as  an  explanation  of  why  sunspots  are  dark. 
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Variations  of  the  Sun’s  radius 

and  temperature 

due  to  magnetic  buoyancy 


Livingston1  has  recently  measured  a  decrease  in  the  surface 
temperature  of  the  Sun  coincident  with  increased  solar  activity. 
He  interpreted  the  temperature  drop  as  implying  a  correspond- 
ing  reduction  in  luminosity.  I  point  out  here  that  surface  cooling 
could  also  be  due  to  a  radial  expansion  of  the  Sun,  with  no 
attendant  reduction  in  luminosity.  There  is  a  plausible  physical 
mechanism  for  such  an  expansion;  namely,  variations  in 
magnetic  buoyancy  due  to  variations  in  the  magnetic  flux  in  the 
convection  zone  over  the  solar  cycle. 

The  solar  luminosity  L  may  be  expressed  in  terms  of  an 
effective  surface  temperature  T  as  L  *  4irR2crT*<  where  R  is 
the  solar  radius.  For  constant  R}  a  decrease  in  T  implies  a 
decrease  in  L.  Alternatively,  there  can  be  a  cooling  (or  heating) 
of  the  surface  due  to  a  radial  expansion  (or  contraction)  of  the 
Sun  with  constant  L.  For  small  changes  A  T  and  A R  with  L 
constant, 


AR  AT 
R  *  2  T 


(1) 


As  T- 6,000  K,  a  drop  in  surface  temperature  of  1  K  would 
correspond  to  a  relative  expansion  A R/R  —  3  x  10  4. 

The  concept  of  magnetic  buoyancy  was  introduced  by  Parker 
and  Jensen3  and  is  an  important  ingredient  in  the  solar  dynamo. 
Consider  an  isolated  magnetic  flux  tube  of  field  strength  B  in  the 
solar  interior.  The  sum  of  the  gas  pressure  p,  and  the  magnetic 
pressure  pm  =  B2/$ir  inside  the  tube  must  balance  the  external 
gas  pressure  p„  so  p,  +  pm  =  pe,  and  p,  <  pe.  If  the  flux  tube  is  in 
thermal  equilibrium  with  its  surroundings  (Tt  =  7U,  then  the 
density  inside  the  tube  is  lower  than  the  density  of  the  surround¬ 
ings  (Pi<p«)  and  the  flux  tube  is  buoyant.  Jensen3  pointed  out 
that,  because  of  the  density  depression  inside  a  magnetic  flux 
tube,  the  presence  of  magnetic  flux  tubes  in  the  solar  interior  will 
cause  the  volume  of  the  Sun  to  increase;  thus,  we  would  expect 
the  radius  of  the  Sun  to  vary  with  the  solar  cycle,  with  maximum 
radius  occurring  near  sunspot  maximum.  Jensen  s  suggestion, 
which  has  attracted  little  previous  attention,  is  pursued  here 

How  much  expansion  and  contraction  of  the  Sun  might  we 
expect  due  to  variations  in  magnetic  buoyancy  over  the  solar 
cycle?  A  precise  estimate  would  require  detailed  knowledge  of 
the  magnetic  field  distribution  in  the  convection  zone  over  the 
solar  cycle  and  a  complete  calculation  of  the  effect  of  the 
magnetic  field  on  the  structure  of  the  solar  envelope.  This  is 
beyond  present  capabilities,  so  instead  a  crude  calculation  is 
used  to  get  a  very  rough  estimate.  Consider  a  toroidal  flux  tube 
that  circles  the  Sun  at  spherical  radius  R,  and  constant  latitude 
4>.  Let  the  tube  have  field  strength  B  and  cross-sectional  radius 
a.  The  mass  of  gas  within  the  tube  is  less  than  would  occupy  the 
same  volume  in  the  absence  of  the  magnetic  field.  The  difference 
in  densities  outside  and  inside  the  tube  is  given  by 
Ap  =pe-p.  =  {Pt-Pi)/RTt  =  pJRTx 


where  T\  *  Tt  =  Tx  is  the  temperature  of  the  flux  tube  and  its 
immediate  surroundings.  The  total  mass  defect  in  the  tube  is 
equal  to  the  volume  of  the  tube  times  Ap,  or 


AM  =  27rRt  cos  <t>  7ra2pmiRTi 

Assuming  that  this  mass  AM  is  spread  over  a  thin  spherical  shell 
of  mean  radius  Ri  and  thickness  A R,  then  the  mass  in  this  shell  is 
given  by 

AM  =  47rR?ARpe 

Equating  the  two  expressions  for  AM  gives 

^  (2) 

R  2  Vr,/\R/  pt 

This  expression  gives  an  estimate  of  the  relative  change  in  radius 
due  to  a  single  toroidal  flux  tube. 


To  obtain  a  numerical  estimate  we  can,  for  example,  use  the 
results  of  the  study  by  Weiss4  of  magnetic  flux  tubes  in  the  solar 
convection  zone,  based  on  equipartition  of  magnetic  energy  and 
kinetic  energy  of  convection.  Weiss  finds  that  a  typical  flux  tube 
has  a  total  flux  of  4  x  1021  Mx.  The  radius  and  field  strength  of 
the  flux  tube  vary  slowly  from  5,000  km  and  5,700  G  at  the 
bottom  of  the  convection  zone  to  14,000  km  and  600  G  at 
photosphere.  The  ratio  pm/p.  varies  greatly,  from  6  x  10~*  at  the 
bottom  of  the  convection  zone  to  2.8  x  10" 1  at  the  photosphere. 
Thus,  the  dominant  contribution  to  A R  will  come  from  flux 
tubes  near  the  top  of  the  convection  zone.  Using  values  given  by 
Weiss  for  a  and  pm/pe  in  equation  (2),  for  R\/R  in  the  range 
0.98-1 .0,  we  obtain  values  of  A  R/R  in  the  range  5  x  10""  cos  <f> 
to  2  x  10  4  cos  0.  For  example,  a  single  flux  tube  at  a  depth  of 
1,000km  at  middle  latitude  gives  AR/R— 2x10  \  Now,  to 
estimate  the  change  in  solar  radius  with  the  solar  cycle,  we  need 
to  estimate  the  difference  in  total  magnetic  flux  in  the  upper 
convection  zone  between  solar  maximum  and  minimum.  If  we 
assume  that  the  total  magnetic  flux  emerging  at  the  solar  surface 
is  a  good  indicator  of  the  total  flux  in  the  uppermost  layers  of  the 
convection  zone,  then  observations  suggest  a  difference  in  total 
flux  of  the  order  1023  Mx  between  solar  maximum  and  mini¬ 
mum.  This  is  equivalent  to  some  25  or  more  of  the  individual  flux 
tubes  discussed  above  and  could  imply  a  relative  change  in 
radius  A  R/R  —  5  x  10  4  or  more  and  a  corresponding  drop  in 
surface  temperature  of  — 1.5  K  or  more. 

Magnetic  field  strengths  in  the  convection  zone  may  be  higher 
than  the  equipartition  limit  used  in  the  estimate  above.  Indeed, 
observations5  have  shown  that  almost  all  of  the  magnetic  flux  in 
the  quiet  photosphere  is  concentrated  to  strengths  of  1,000- 
2,000  G,  well  above  the  equipartition  limit,  and  the  theory  of 
Galloway,  Proctor  and  Weiss6  predicts  that  magnetic  fields  will 
be  concentrated  to  10  times  the  equipartition  limit  or  more  by 
the  convection.  Equation  (2)  can  be  rewritten  as 

A  R  ^  FB  cos  4> 

R  \6rrRR\p e 

where  F  =  ira2B  is  the  total  magnetic  flux  in  the  tube.  For  a  tube 
of  fixed  total  flux,  the  relative  change  in  radius  is  proportional  to 
the  magnetic  field  strength.  Thus,  higher  field  strengths  in  the 
convection  zone  would  increase  the  estimate  of  A  R/R  for  a  fixed 
difference  in  total  flux  over  the  solar  cycle.  With  higher  field 
strengths,  flux  tubes  deeper  in  the  convection  zone  would 
contribute  significantly  to  A  R/R.  Although  there  is  uncertainty 
about  magnetic  fields  in  the  convection  zone,  the  estimates  do 
suggest  that  the  surface  cooling  of  -6  K  observed  by  Livingston 
may  be  at  least  partly  caused  by  expansion  due  to  magnetic 
buoyancy. 

The  historical  record  of  measurements  of  the  solar  radius 
gives  no  clear  or  consistent  picture7  9 .  The  currently  accepted 
value10  (R  =  959.63"  at  1  ad)  dates  back  to  1891  and  is  subject 
to  correction  for  irradiation11.  Variability  of  R  with  amplitudes 
ranging  from  —0.05”  to  —2"  and  periods  of  7, 8, 1 1  and  22  yr  has 
been  reported.  Giannuzzi12,  for  example,  found  variations  in  R 
with  amplitude  0.2-0. 5"  and  period  22  yr.  Meyermann13  repor¬ 
ted  variations  with  total  amplitude  -0.15"  (A R/R  -  1.5  x  10  4) 
and  period  1 1  yr  with  maximum  radius  occurring  very  near 
sunspot  maximum;  this  is  in  good  agreement  with  the  theoretical 
argument  above.  Gething*  finds  variations  in  R  with  amplitude 
—0.1 -0.2",  but  with  no  obvious  correlation  with  the  solar  cycle. 
Improved  measurements  of  the  solar  radius  over  a  solar  cycle 
are  needed  to  test  the  ideas  expressed  here.  Recent  improve¬ 
ments9 14  in  the  theoretical  definition  of  the  solar  radius  should 
give  a  considerable  increase  in  accuracy.  Changes  in  opacity 
during  expansion  may  affect  the  visual  change  in  solar  radius. 
Note  that  the  solar  oblateness  measured  by  Dicke  and  Golden- 
berg15  (A R/R  -  5  x  10  '*)  and  later  disputed  by  Hill  and  Steb- 
bins16  is  an  order  of  magnitude  smaller  than  the  variation  in 
mean  radius  estimated  above.  Even  if  the  expansion  due  to 
magnetic  buoyancy  is  too  small  to  produce  a  measurable  change 
in  surface  temperature,  it  may  still  be  large  enough  to  show  up 
directly  as  a  change  in  radius.  Accurate  monitoring  of  the  solar 


radius  over  the  solar  cycle  would  give  a  measure  of  the  variation 
of  total  magnetic  flux  in  the  convection  zone  and  add  to  our 
understanding  of  the  solar  dynamo.  Measurements  of  the  solar 
radius  from  space  are  particularly  desirable. 
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Abstract 

Umbral  oscillations  in  sunspots  are  identified  as  a 
resonant  response  of  the  umbral  atmosphere  to  forcing  by 
oscillatory  convection  in  the  subphotosphere.  The  full, 
linearized  equations  for  magneto-atmospheric  waves  are  solved 
numerically  for  a  detailed  model  of  the  umbral  atmosphere,  for 
both  forced  and  free  oscillations.  Resonant  ,rfast"  modes  are 
found,  the  lowest  mode  having  a  period  of  153  s,  typical  of 
umbral  oscillations.  A  comparison  is  made  with  a  similar 
analysis  by  Uchida  and  Sakurai  (1975)  ,  who  calculated  resonant 
modes  using  an  approximate  ( "quasi-Alfven" )  form  of  the  wave 
equations.  Whereas  both  analyses  give  an  appropriate  value 
for  the  period  of  oscillation,  several  new  features  of  the 
motion  follow  from  the  full  equations. 
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1.  Introduction 


Observations  of  umbral  oscillations  in  sunspots  (Beckers 
and  Schultz  1972;  Bhatnagar  and  Tanaka  1972;  Giovanelli  1972; 
Rice  and  Gaizauskas  1973;  Phillis  1975;  Moore  and  Tang  1975) 
give  a  fairly  consistent  picture  of  the  oscillations  as  a 
resonant  wave  mode  in  the  umbra.  Doppler  shifts  in  Ha  and  in 
photospheric  lines  show  a  continuous  periodic  oscillation  in 
sunspot  umbrae,  with  a  fairly  well  defined  period  in  the  range 
145-190  s. 

Moore  (1973)  has  argued  that  the  driving  mechanism  for 
umbral  oscillations  is  overstable  oscillatory  convection  in  a 
shallow  subphotospheric  layer  in  the  umbra  (see  also  Mullan  and 
Yun  1973) .  Uchida  and  Sakurai  (1975)  also  consider  the  driving 
mechanism  to  be  overstable  convection,  and  they  interpret  the 
umbral  oscillations  as  a  resonant  response  of  the  umbral 
atmosphere  to  this  forcing.  Their  resonant  mode  is  a  standing 
quasi-Alfven  wave  trapped  in  the  photosphere  and  chromosphere. 
In  a  similar  approach,  Antia  and  Chitre  (1979)  identify  the 
umbral  oscillations  with  an  overstable  fast  magneto-atmospheric 
wave  mode  in  the  umbra.  Their  calculations  account  more  fully 
for  the  effects  of  compressibility. 

In  this  paper,  we  adopt  a  point  of  view  very  close  to  that 
of  Uchida  and  Sakurai  (1975).  We  shall  argue  that  the  umbral 
oscillations  are  a  resonant  response  to  overstable  convection 
in  a  thin  subphotospheric  layer  and  that  the  resonant  wave  mode 


is  a  fast  magneto-atmospheric  wave.  Our  basic  model  atmosphere 
is  similar  to  Uchida  and  Sakurai's.  There  are,  however,  several 
fundamental  differences  between  our  work  and  Uchida  and  Sakurai' 
which  lead  to  a  new  understanding  of  the  umbral  oscillations. 
These  differences  include  the  following: 

(1)  We  do  not  adopt  the  approximation  that  leads  to  the 
quasi-Alfven  wave;  rather,  we  solve  the  complete  linearized 
magneto-atmospheric  wave  equations.  We  show  that,  although 
the  quasi-Alfven  approximation  yields  a  good  estimate  of  the 
period  of  oscillation,  it  fails  to  describe  certain  important 
features  of  the  motion. 

(2)  Instead  of  assuming  a  reflecting  lower  boundary  at 
the  base  of  the  photosphere,  we  include  a  semi-infinite  lower 
layer  in  our  model  atmosphere  to  represent  the  umbral  convection 
zone.  We  show  that,  as  far  as  umbral  oscillations  are  concerned 
this  lower  layer  acts  much  like  a  reflecting  boundary,  but  not 
for  the  reasons  given  by  Uchida  and  Sakurai. 

(3)  In  addition  to  calculating  free  eigenmodes  of  oscill¬ 
ation  in  our  model  atmosphere,  we  also  calculate  the  response 
to  forcing  at  different  frequencies. 

(4)  We  show  that  the  downward  reflection  of  wave  energy 
is  not  total,  but  that  a  small  fraction  of  energy  escapes  into 
the  corona  in  the  form  of  acoustic  waves  along  the  magnetic 
field  lines. 


2.  Basic  Equations  and  the  Uinbral  Model 


In  our  simplified  treatment  of  a  sunspot  umbra,  we  consider 
the  undisturbed  magnetic  field  to  be  uniform  and  vertical, 

B  =  (0,0, Bq) .  The  atmosphere  is  assumed  to  be  a  compressible, 
inviscid,  perfectly  conducting  gas  under  uniform  gravity 
g(=  0.274  km  s  )  in  the  negative  z-direction.  The  undisturbed 
pressure  p(z),  density  p(z),  and  temperature  T(z)  are 
functions  of  height  z  only,  and  hydrostatic  equilibrium 
requires  that 
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We  then  consider  small  adiabatic  perturbations  of  this  equilibrium 
atmosphere.  There  is  no  preferred  horizontal  direction,  so  in 
cartesian  coordinates  we  may  take  the  horizontal  wavenumber  in 
the  x-direction  and  assume  that  the  perturbation  velocity 
u  =  (u,v,w)  has  the  form  u  =  u  (z)  exp  [i  (kx  -  oat)  ]  ,  with 
u(z)  =  [u(z),  v(z),  w(z)].  Starting  with  the  linearized  equations 
of  induction  and  conservation  of  mass,  momentum,  and  energy,  we  can 
eliminate  the  perturbations  in  pressure,  density,  and  magnetic 
field  and  arrive  at  a  single  vector  equation  for  the  perturbation 
velocity  (cf .  Ferraro  and  Plumpton  1958)  .  The  y  component  of 
this  equation  is 


r  2  d2  .  2.-  A 

[vA  — ^  +  0)  ]v  =  0  , 

dz 


(1) 


2  *5 

where  =  (Bo/4np (z) )  is  the  Alfven  speed.  This  equation 

is  decoupled  from  the  other  two  components  and  represents  a 
pure,  transverse  Alfven  wave  with  horizontal,  incompressible 
motion.  The  x  and  z  component  equations  are 

t v2(-^-  -  k2)  -  c2k2  +  to2] G  +  iktc2  -  g]w  =  0,  (2) 

dz 

2 

[C2  -^2  -  yg  ^  +  u>2lw  +  ik[c2  ^  -  (y-l)glu  =  0,  (3) 

dz 

where  c  =  (yRT  (z)  )**  is  the  sound  speed  and  y  is  the  ratio 
of  specific  heats.  Equations  (2)  and  (3)  represent  fully 
coupled  magneto-atmospheric  waves  in  which  compression,  buoyancy, 
and  magnetic  forces  all  play  a  role. 

Equations  (2)  and  (3)  also  describe  the  vertical  dependence 
of  waves  in  cylindrical  coordinates  (r,0,z).  If  we  take  the 
perturbation  velocity  u  *  (u  ,ufl,u  )  in  the  form 

*w  IT  V  Z 

v»r  =  u(z)  [kJm+1(kr)  -  ™  Jm(kr)  ]exp[i(m9  -  tot)],  (4) 

J  (kr) 

Uq  =  -im  u(z)  -  exp [i  (m0  -  cot)  ]  ,  (5) 

u  =  -ik  w (z)  J  (kr)  exp  [i  (m0  -  cot )  ]  ,  (6) 

z  m 

we  arrive  at  equations  identical  to  (2)  and  (3) ,  except  that  in 
this  case  k  is  a  radial  wavenumber.  We  shall  consider  only 
axisymmetric  modes  (m=0) ,  for  which  Ug  =  0  and 


(7) 


ur  =  ku (z) (kr) exp (-iwt) , 

u  =  -ik  w(z) J  (kr)exp(-iwt) .  (8) 

(There  are  also  axisymmetric  modes  with  u0  ?  0,  representing 

pure,  torsional  Alfven  waves;  these  are  not  represented  by 

the  form  (4) -(6).)  If  we  require,  as  Uchida  and  Sakurai  do, 

that  the  radial  velocity  ur  vanish  at  the  edge  of  the  umbra, 

r  =  a,  then  there  is  a  discrete  set  of  radial  wavenumbers  k^ 

defined  by  (k ^a)  =0.  We  shall  consider  only  the  first  zero 

of  J^,  ka,=  3.83,  and  take  the  same  value  k  ^  =  1100  km  as 

Uchida  and  Sakurai,  corresponding  to  an  umbral  radius 
■% 

a  =  4200  km. 

The  temperature  distribution  T(z),  and  hence  the  density 
distribution  p(z),  must  be  specified  to  represent  the  umbral 
atmosphere.  We  adopt  as  our  model  umbra  the  three-layer 
atmosphere  shown  in  Figure  1.  This  same  atmosphere  was  used  in 
a  previous  study  of  reflection  of  Alfven  waves  in  the  umbra 
(Thomas  1978) ,  and  the  upper  two  layers  are  identical  to  the 
model  atmosphere  of  Uchida  and  Sakurai.  Layer  2,  representing 
the  umbral  photosphere  and  chromosphere,  is  isothermal  at 
temperature  Tq.  Layer  3  represents  the  corona  and  is  isothermal 
at  temperature  T^  (>Tq)  •  The  transition  region  is  represented 
as  a  discontinuity  in  temperature  (and  density)  at  height  z  =  zfc.' 
Layer  1,  which  represents  the  umbral  convection  zone,  is  assumed 
to  have  a  linear  temperature  distribution  T(z)  =T  -  0z 
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with  8>0.  The  value  of  the  temperature  gradient  8  may 

be  chosen  equal  to  or  slightly  greater  than  the  adiabatic 

gradient  8  =  g/c  .  Let  a  =  H  8/T_  be  a  nondimensional 
s  p  o  o 

measure  of  the  temperature  gradient  8  in  layer  1,  where 
Hq =  RTQ/g  is  the  density  scale  height  in  layer  2.  The 
temperature,  density,  sound  speed,  and  Alfven  speed  in  each 
layer  are  then  given  as  follows. 

Layer  1  (z  <  0) : 


=  Tc(l  -  2£)  ,  H0-BT0/9. 

o 


p(z)  =  Po(l  -  ^)  a  , 
o 

(9) 

2  .  x  2  . .  az  x  2 

C  (Z)  =  C0(l  -  |p)  #  cq 

o 

=  yhto, 

-  n  Ot-1 

VA<2>  =  VAo(1-^  “  ' 

v2  =  B2/4irp  . 

Ao  o  o 

Layer  2  (0  <  z  <  zfc)  : 


T(z)  =  To,  Ho  =  RTo/g, 

P (z)  =  P  exp ( -z/H  ) , 

°  °  (10) 
c2(z)  =  c2  =  YRTq, 

VA(z)  =  VAo  exp(z/«0),  v2q  =  B2/4ttpo. 
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Layer  3  (z  >  zfc) : 


T  (z) 


=  Hx  =  RTj/g 


c2  (z) 


2,  x 

VA  2 


Po(T^)expt_Zt/Ho  ~  (z-z^/H^, 


C1  =  YRTr 
T 

vAo(r)wlzt/Ho +  <*-*t>/V  • 
o 


We  note  here  that,  although  we  do  compute  solutions  of 
the  wave  equations  (2)  and  (3)  in  the  entire  three-layer 
atmosphere,  we  also  present  results  for  simpler  versions  of 
the  atmosphere,  where  the  lower  layer  is  replaced  by  a  reflect¬ 
ing  boundary  or  where  the  upper  layer  is  omitted  and  layer  2 
extends  to  z  -  00 . 

For  comparison  with  observations,  we  adopt  the  following 

numerical  values  for  the  parameters  in  our  umbral  model: 

Tq  =  4500  K,  =  2X106  K,  Bq =  1000  G,  Pq =  5*10-7  g  cm'3, 

y  =  5/3,  and  a  =  0.5.  The  sound  speed  and  density  scale  height 

in  layer  2  are  then  c  =7.9  km  s  1  and  H  =136.5  km,  and 

o  o 

the  Alfven  speed  at  z  =  0  is  vAq  =  4.0  km  •  We  take  the 
height  of  the  transition  region  to  be  =  20  Hq=  2730  km. 


3.  Analysis 


In  order  to  compute  wave  modes  in  our  model  umbral  atmos¬ 
phere,  we  need  to  solve  the  basic  coupled  wave  equations  (2) 
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and  (3)  for  two  specific  cases,  uniform  temperature  (layers  2 

and  3)  and  linearly  varying  temperature  (layer  1) .  Consider 

first  the  isothermal  case.  The  problem  of  magneto-atmospheric 

waves  in  an  isothermal  atmosphere  with  a  uniform  vertical 

magnetic  field  was  first  treated  by  Ferraro  and  Plumpton  (1958) 

in  an  important  early  paper  that  has  escaped  notice  in  subsequent 

treatments  of  this  problem  by  Uchida  and  Sakurai  (1975)  and 

Hollweg  (1979) .  Ferraro  and  Plumpton  solve  the  full  wave 

equations  (2)  and  (3) ,  whereas  Uchida  and  Sakurai  and  Hollweg 

solve  approximate  forms  of  the  equations  which  hold  only  in 
2  2 

the  case  c  <<  v^.  This  approximation  is  not  good  m  the  low 
photosphere,  so  we  prefer  to  use  the  full  equations.  In 
contrast  to  the  case  of  a  uniform  horizontal  magnetic  field 
(Nye  and  Thomas  1976) ,  the  wave  equations  for  an  isothermal 
atmosphere  with  a  uniform  vertical  magnetic  field  do  not  have 
solutions  expressible  in  terms  of  tabled  special  functions. 
Ferraro  and  Plumpton  construct  power  series  solutions  to  the 
equations  and  give  only  limited  numerical  results.  For  our 
purposes,  it  is  more  convenient  to  rely  on  direct  numerical 
integration  of  the  equations. 

It  is  useful  to  have  asymptotic  solutions  of  the  wave 
equations  for  large  z.  Suppose  the  transition  zone  is  absent 


L 


2  2 

and  layer  2  extends  to  z-°°.  Then,  for  large  z,  v^  >>  c  , 


and  the  asymptotic  solutions  of  Equations  (2)  and  (3)  for 
z ->  oo  are  (cf.  Ferraro  and  Plumpton  1958) 
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-■  --  -  —  - - 


J 


f  ! 


u(z)  -  A  exp(-Kz/HQ)  +B  exp  (-z/2Hq)  exp  [±i  (ft^-1*)  ^z/H^  ,  (12) 


w(z)  -  C  exp  (z/2Hq)  exp  [±i  (ft2-1*)  ^z/H^]  , 


(13) 


where  K  =  H  k  is  a  nondimensional  horizontal  wavenumber  and 
o 

ft  =  Hq<jo/co  is  a.  nondimensional  frequency.  For  ft  >  *5  the  plus 
sign  in  the  exponents  corresponds  to  an  upward  propagating 
wave  and  the  negative  sign  to  a  downward  propagating  wave. 

The  acoustic  cutoff  frequency  corresponds  to  ft  = H ,  and  for 
ft  <  *3  the  asymptotic  solution  represents  purely  evanescent 
behavior. 

The  case  of  a  linear  temperature  variation  (layer  1)  has 
% 

not  been  studied  previously.  We  have  not  found  a  useful 
representation  of  the  solution  to  the  wave  equations  in  terms 
of  special  functions,  so  we  again  use  direct  numerical 
integration . 

In  analyzing  our  computed  solutions,  it  is  useful  to 
consider  the  distribution  of  wave  energy  density  with  height. 
The  total  energy  density  E  is  given  by 


E  =  *ipu2  +  »jpw2  +  -^-[K2u2  +  (w-H  2  +  (l^i)  2w2] 

ft*  o  az  y 

2  2 

,  i^/  °  \^2  .  .  Ar„2/dUv  2  .  T,2~2. 

+  hpl— jJw  +  >,P  -ylH  (gj)  +K  u  ] 

ygft  c 


(14) 


in  layer  2, 


N 


2 


'  Y  ,H1 


(Y-l) 


in  layer  3 ,  and 


N2  = 


Y-gy-l 


Yd  -*r>J 


_g_ 

H 


in  layer  1.  The  five  terms  on  the  right  hand  side  of  equation 
(14)  represent ,  in  order,  the  kinetic  energy  density  due  to 
horizontal  motions,  the  kinetic  energy  density  due  to  vertical 
motions,  the  ^potential  energy  density  associated  with  adiabatic 
compression  or  expansion,  the  potential  energy  density  associated 
with  the  buoyant  force,  and  the  potential  energy  density  stored 
in  the  magnetic  field  perturbation. 

The  wave  equations  (2)  and  (3)  are  solved  numerically 
using  a  generalized  Newton-Raphson  relaxation  method.  The 
equations  are  written  as  a  set  of  linear,  first-order  ordinary 
differential  equations  with  the  velocities  (and  eigenvalues, 
if  desired)  expressed  as  some  trial  value  plus  a  correction. 

They  are  then  linearized  with  respect  to  the  corrections.  The 
atmosphere  is  divided  into  a  grid  of  equally  spaced  points. 

Grid  spacings  between  0.1  Hq  and  0.25  Hq  give  sufficient 
accuracy  in  all  of  the  calculations.  Using  a  finite  difference 
scheme,  the  calculation  is  started  at  one  boundary,  and,  using 
the  equations,  the  corrections  to  the  initial  trial  values  are 
calculated  at  each  internal  grid  point.  The  conditions  at  the 


other  boundary  close  the  set  of  equations  for  the  necessary 
corrections  at  every  mesh  point.  The  corrections  are  then 
added  to  the  velocities  (and  eigenvalues)  giving  the  trial 
quantities  for  the  next  iteration.  The  process  is  repeated 
until  the  maximum  correction  divided  by  the  associated  quantity 
is  less  than  some  pre-assigned  value  (0.1%). 

4 .  Forced  Oscillations 

To  simulate  the  forcing  by  overstable  convection  in  a 
thin  subphotospheric  layer  we  apply  an  oscillation  of  fixed 
frequency  and  horizontal  wavenumber  and  of  unit  amplitude  on 
the  plane  z  =  0.  We  consider  separately  the  forcing  due  to 
a  purely  horizontal  motion  (u(0)  -1,  w(0)  =0)  and  a  purely 
vertical  motion  (G(0)  =0,  &(0)  =1)  at  z  =  0.  A  true  repre¬ 
sentation  of  the  forced  oscillations  will  then  be  some  linear 
combination  of  these  two  forced  solutions. 

We  shall  consider  these  forced  solutions  in  two  stages. 
First  we  consider  only  a  two-layer  atmosphere  consisting  of 
layers  1  and  2,  with  layer  2  extending  to  z  = 00  and  with 
layer  3  (the  corona)  absent  (i.e.,  z  -*•<»)  .  Then  we  consider 
solutions  for  the  full  three-layer  atmosphere.  This  approach 
will  elucidate  the  role  of  the  chromosphere  -  corona  transition 
in  providing  downward  reflection  of  wave  energy. 

4.1  TWO-LAYER  MODEL 

Here  layer  2  is  semi-infinite  (z .•♦«>)  and  layer  3  (the 


corona)  is  absent.  In  layer  2  we  assume  that  the  asymptotic  solu 

tion  (12),  (13)  is  valid  above  z  =  20  H  and  use  this  as  an 

o 

end  point  of  our  numerical  integration.  At  z  ■  20  Hq  we 
apply  an  outward  radiation  condition ,  that  is,  we  allow  only 
outward  propagation  (in  the  positive  z-direction)  of  energy. 

This  leads  to  the  following  matching  conditions  at  z  =  20  Hq: 


du 

dz 


exp(-Kz/Ho)  +  -*5]exp(-z/2Ho)exp[i(n2-»,)l5z/Ho] 


2  L  1 

exp(-z/2H  )exp[i(£2  -1*)  z/H  ]  exp(-Kz/H  ) 

O  O  K  O 


u  =  0, 


(15) 


-  -pma2-'*)15  +  jj]w  =  o. 

dz  H 

o 


The  solution  in  layer  1  is  more  difficult  to  obtain  because 
there  is  no  simple  asymptotic  form  in  the  region  of  interest. 

We  obtain  a  simple  condition  of  outward  propagation  (in  the 
negative  z-direction)  by  considering  the  kinetic  energy 
density,  which  must  remain  finite  as  z  -*■  Since  the 

density  is  increasing  linearly  (for  a =  h) ,  the  velocity 
components  must  drop  off  at  least  as  fast  as  z  ** .  Consequently, 
the  derivatives  must  drop  off  even  faster  so  we  apply  the 
conditions 


i 


du 

dz 


0  , 


dw 

dz 


0 


(16) 


at  some  depth  z  =  -h.  To  insure  that  the  error  introduced  by 
this  is  smaller  than  our  allowed  numerical  error  (0.1%),  the 
value  of  h  is  increased  until  changes  in  h  have  no  effect 
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on  u  and  w.  We  found  such  a  suitable  value  to  be  h = -20  Hq. 

In  summary,  we  evaluate  numerically  the  solution  for 
either  horizontal  (u(0)  =1,  w(0)  =0)  or  vertical  (u(0)  =0, 
w(0)  =1)  forcing,  subject  to  the  outward  radiation  conditions 
(15)  and  (16),  applied  at  z =  20  Hq  and  z = -20  Hq, 
respectively.  In  the  resulting  solutions,  du/dz  and  dw/dz, 
and  hence  the  pressure  and  magnetic  field  perturbations,  are 
discontinuous  across  z  =  0.  This  is  acceptable  in  view  of 
the  fact  that  the  plane  z  =  0  represents  an  overstable  layer 
of  finite  thickness. 

In  the  case  of  horizontal  forcing  (u(0)  =1,  w(0)  =0), 
we  have  computed  solutions  for  a  fixed  wavenumber  K= 0.124 
(corresponding  to  k  ^  =  1100  km  and  Hq  =  136.5  km)  and  for 
many  different  frequencies.  Figure  2  shows  a  plot  of  the  peak 
energy  density  (which  occurs  at  some  z  s  0)  of  the  forced 
oscillation  as  a  function  of  the  forcing  frequency  ft.  This 
graph  clearly  shows  the  existence  of  a  resonant  response  at 
certain  frequencies.  Numerical  values  of  the  first  three 
resonant  frequencies  ftR  determined  in  Figure  2  are  listed 
in  Table  1,  along  with  the  corresponding  dimensional  periods 
using  the  numerical  values  of  atmospheric  parameters  listed 
in  Section  2.  Each  of  these  frequences  lies  above  the  acoustic 
cutoff  frequency,  ftc  =  1/2  (ojc  =  cq/2Ho)  and  the  resonant 
modes  are  "fast"  magneto-atmospheric  waves.  Similar  calcula¬ 
tions  for  the  case  of  vertical  forcing  show  a  resonant  response 


at  the  same  resonant  frequencies  as  in  the  case  of  horizontal 
forcing.  The  period  of  the  lowest  resonance,  153  s,  lies 
comfortably  in  the  range  of  observed  periods  of  umbral 
oscillations . 

Figure  3  shows  the  solution  for  the  horizontal  velocity 
u(z)  as  a  function  of  height  in  layer  2  for  each  of  the  first 
three  resonant  frequencies.  The  computed  magnitude  is  finite 
in  each  case  because  the  resonant  frequency  is  only  approxi¬ 
mately  matched.  These  first  three  resonant  modes  closely 
resemble  the  corresponding  free  eigenmodes  computed  by 
Uchida  and  Sakurai  (1975,  Figure  3)  for  their  model;  we 
shall  say  more  about  this  correspondence  in  Section  5.  Each 
higher  resonance  has  an  additional  node  in  the  solution  for 
u(z) .  The  numerical  solutions  for  u(z)  agree  with  the 
asymptotic  solution  to  within  0.1%  for  z  > 6Hq;  this  more 
than  justifies  our  matching  to  the  asymptotic  solution  at 

z  =  20H  . 
o 

The  horizontal  velocity  u(z)  tends  toward  zero  as 
z  -*■<»,  due  primarily  to  the  exponentially  increasing  Alfven 
speed.  The  vertical  velocity,  however,  grows  with  height, 
eventually  increasing  exponentially  according  to  its 
asymptotic  form  (13) .  Because  of  this  rapid  growth  it  is 
more  convenient  to  plot  the  logarithmic  derivative 
HQ(dw/dz)/w  than  w  in  this  layer.  Figure  4  shows  the 
computed  magnitude  of  this  logarithmic  derivative  as  a  function 
of  height  is  the  upper  layer,  for  horizontal  forcing  at  the 


lowest  resonant  frequency.  We  see  that  | H^tdw/dz) /w| 

attains  its  asymptotic  value  ft  (cf.  Equation  (13))  for 

z  *  10H  . 
o 

The  behavior  of  the  vertical  velocity  in  this  case  of 
purely  horizontal  forcing  is  of  some  interest.  Although 
vertical  motions  are  not  forced  directly  (0(0)  =  0) ,  they 
are  produced  above  and  below  z  =  0  by  the  "squeezing"  effect 
of  the  divergent  horizontal  motions.  In  the  upper,  isothermal 
layer  vertical  motions  can  propagate  upward  in  the  form  of  a 
pure  acoustic  wave  along  the  magnetic  field  lines.  Since  the 
sound  speed  is  constant  in  the  upper  layer,  there  is  no 
downward  reflection  of  such  an  acoustic  wave,  and  it  will 
propagate  indefinitely  upward  with  growing  amplitude  due  to 
decreasing  density.  This  provides  a  mechanism  for  escape  of 
energy,  which  means  we  are  not  dealing  with  a  perfect  resonator. 
For  horizontal  forcing  at  our  chosen  frequency  and  wavenumber, 
however,  only  a  small  fraction  of  the  energy  is  converted  into 
acoustic  energy  that  escapes  to  z  =  °°.  Figure  5  shows  the 
distribution  of  total  energy  density  (Equation  (14) )  with 
height  for  the  lowest  resonance.  The  normalized  curve  is 
identical  for  horizontal  and  vertical  forcing.  For  positive 
z  the  energy  density  drops  sharply  at  first  and  then  levels 
off  rather  abruptly  at  z  ~  5Ho  to  a  nearly  constant  value  of 
about  2%  of  the  peak  value.  By  considering  individual  terms 
in  the  expression  for  total  energy.  Equation  (14),  we  have 
established  that  nearly  all  the  energy  near  z  =  0  is  in  the 


form  of  magnetic  energy  and  kinetic  energy  due  to  horizontal 


motion,  whereas  for  z  ;>  5Hq  the  energy  is  nearly  all  due  to 
vertical  motions,  buoyancy,  and  compression.  The  sharp  drop 
in  energy  density  between  z  =  0  and  z  =  5Hq  is  due  to 
strong  downward  reflection  by  the  exponentially  growing  Alfven 
speed.  The  small,  constant  energy  density  above  z  =  10Hq 
represents  escaping  energy  in  the  form  of  a  pure  acoustic 
wave  along  the  vertical  magnetic  field  lines. 

For  forcing  at  the  resonant  frequency,  the  velocities 
and  energy  density  in  layer  1  (z  <  0)  are  much  smaller  than 
in  layer  2  (z  >  0)  and  are  therefore  shown  separately. 

Figure  6  shows  the  horizontal  and  vertical  velocities  and 
the  total  energy  density  in  the  lower  layer  for  both  horizontal 
and  vertical  forcing  at  the  lowest  resonant  frequency.  For 
horizontal  forcing,  the  horizontal  velocity  drops  off  rapidly 
with  depth  and  the  vertical  velocity  rises  to  a  maximum  at 
z  ~ -2Ho  before  decaying  rapidly  as  z-*--00.  The  purely 
horizontal,  compressive  (kj*0)  motion  at  z  =  0  induces 
vertical  motions  for  z  <  0  through  a  "squeezing"  action. 

The  resulting  compressive  vertical  motions  below  z  =  0  have 
an  acoustic  character  and  are  reflected  strongly  upward  by 
the  increasing  sound  speed.  For  vertical  forcing,  horizontal 
motions  are  produced  below  z  =  0  by  a  similar  "squeezing" 
effect,  but  to  a  lesser  extent.  The  total  energy  in  the  lower 
layer  is  less  for  the  case  of  vertical  forcing. 

The  buoyant  force  is  not  important  in  causing  upward 


6 


reflection  in  layer  1.  Solutions  for  the  case  a  =  (Y_l)/Y  =  0.4, 
corresponding  to  an  adiabatic  temperature  gradient  (neutral 
stability)  in  layer  1,  have  the  same  character  as  the  solution 
in  Figure  6,  for  which  a =0.5  (slightly  superadiabatic) . 

For  forcing  at  the  lowest  resonant  frequency,  the  total 
energy  density  in  the  lower  layer  is  six  orders  of  magnitude 
smaller  than  in  layer  2,  and  hence  doesn't  even  appear  in 
Figure  5.  The  resonant  modes  are  almost  totally  confined  to 
layer  2.  To  illustrate  the  sharpness  of  the  resonant  response, 
it  is  useful  to  look  at  the  response  to  forcing  at  a  frequency 
slightly  off  resonance.  Figure  7  shows  the  distribution  of 
total  energy  density  for  both  horizontal  and  vertical  forcing 
at  frequency  ft=  .743,  about  4%  above  the  lowest  resonant 
frequency.  For  horizontal  forcing,  the  normalized  energy 
density  distribution  in  the  upper  layer  is  nearly  identical 
to  that  of  the  resonant  response  (Figure  5) ,  but  the  actual 
numerical  values  are  six  orders  of  magnitude  smaller  and 
there  is  comparable  total  energy  in  layer  1  (z  <  0) .  For 
vertical  forcing  at  this  off-resonance  frequency,  the  energy 
density  only  drops  to  36%  of  its  peak  value  as  z-*°° ,  and 
there  is  considerably  more  outward  leakage  of  energy  than  in 
the  resonant  response.  In  this  case,  we  expect  the  presence 
of  the  corona  (layer  3)  to  have  an  important  influence.  For 
z>  5Hq,  the  energy  density  is  almost  completely  due  to 
acoustic-like  motions  along  the  magnetic  field  lines;  such 
motions  will  be  reflected  strongly  at  the  chromosphere-corona 
transition . 


62 


4.2  THREE-LAYER  MODEL 

In  Section  4.1  above  we  have  established  that  there  is  a 
resonant  response  of  the  umbral  atmosphere,  at  a  frequency 
appropriate  to  umbral  oscillations,  in  a  two-layer  model 
atmosphere  without  a  corona.  The  necessary  downward  reflection 
is  provided  by  the  rapidly  increasing  Alfven  speed  in  layer  2 
(the  photosphere  and  chromosphere) .  Indeed,  the  energy  density 
of  the  resonant  response  has  dropped  to  2%  of  its  peak  value 
at  the  assumed  height  of  the  transition  region,  z  =  20  Hq 
(Figure  5).  Thus,  for  umbral  oscillations,  we  find  the 
presence  of  the  high-temperature  corona  to  be  of  only  minor 
importance.  This  is  in  sharp  contrast  with  Uchida  and 
Sakurai  (1975)  who  place  heavy  emphasis  on  downward  reflection 
of  energy  at  the  transition  region.* 

For  completeness,  we  consider  the  effect  of  the 
chromosphere-corona  transition  region  on  the  lowest  resonant 
response.  We  have  calculated  the  response  to  horizontal  and 
vertical  forcing  in  the  full  three-layer  model  atmosphere 
(Section  2),  with  the  transition  region  at  z^ = 20  Hg.  The 
solution  in  layer  3  (the  corona)  is  subjected  to  the  outward 
radiation  condition  (15) ,  with  Hq  replaced  by  at 

z  =  40  Hq.  Figure  8  shows  the  distribution  of  total  energy 
density,  which  is  the  same  for  both  horizontal  and  vertical 

*Their  emphasis  on  reflection  at  the  transition  region  is 
unwarranted,  however;  an  investigation  of  their  analytical 
solution  shows  that  almost  all  of  the  reflection  occurs  below 
the  transition  region,  consistent  with  our  model. 


forcing.  A  comparison  of  Figures  5  and  8  shows  that  the 
small  fraction  of  wave  energy  at  z  =  zfc  =  20  HQ  is  almost 
totally  reflected  downward  by  the  sharp  temperature  rise. 

The  energy  density  for  z  >  20  HQ  is  three  orders  of  magni¬ 
tude  smaller  than  at  z <  20  H  and  doesn't  show  up  on  the 

o 

scale  of  these  graphs.  The  downward  reflection  at  z  = 
also  sets  up  a  standing  wave  pattern  that  leads  to  the  small 
oscillations  in  energy  density  with  height  above  z=5  Hq, 
visible  in  Figure  8. 

Downward  reflection  at  the  transition  region  is  of  more 
importance  in  the  case  of  forcing  at  of f -resonance  frequencies, 
where  a  higher  fraction  of  the  peak  wave  energy  penetrates 
to  the  height  of  the  transition. 

5.  Free  Eigenmodes 

The  results  of  our  investigation  of  the  resonant  response 
of  the  umbral  atmosphere  to  forcing  at  z  » 0  allow  us  to  take 
a  simpler  approach  to  the  problem.  We  have  seen  that  the 
forced  resonant  modes  are  strongly  reflected  in  layer  1  (z  <  0) , 
with  negligible  wave  energy  in  that  layer.  Thus,  at  least  as 
far  as  the  resonant  modes  are  concerned,  it  is  reasonable  to 
replace  layer  1  by  a  rigid,  perfectly  conducting  lower  boundary 
at  z=0.  Uchida  and  Sakurai  (1975)  assume  such  a  lower 
boundary  in  their  analysis  of  quasi-Alfven  waves;  however, 
their  stated  justification  for  this  is  quite  different  than 
ours.  Uchida  and  Sakurai  argue  that  the  region  below  z  =  0 
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acts  as  an  almost  rigid  boundary  because  of  the  higher  density 
there.  However,  Thomas  (1978)  has  shown  that  pure  Alfven  waves 
are  only  weakly  reflected  when  propagating  downward  into  a 
region  identical  to  layer  1;  thus,  Uchida  and  Sakurai's 
reasons  for  assuming  a  rigid  lower  boundary  are  questionable. 

In  our  present  calculations  with  the  complete  wave  equations, 
including  the  effects  of  compression,  we  find  that  it  is  the 
compressive  nature  of  the  wave  motions  and  the  increasing 
temperature  (and  sound  speed)  with  depth  that  lead  to  strong 
upward  reflection  from  layer  1. 

Our  results  for  forced  resonances  also  show  that  there 
is  strong  downward  reflection  in  the  isothermal  layer  (layer  2) , 
with  only  a  very  small  fraction  of  the  wave  energy  reaching  as 
high  as  the  chromosphere-corona  transition  region.  This 
suggests  that,  instead  of  an  outward  radiation  condition  at 
some  height,  we  could  simply  assume  another  perfectly  reflect¬ 
ing  boundary  at  some  sufficiently  large  height  z  =  h  without 
changing  the  nature  of  the  wave  mode  appreciably. 

As  a  simplified  model  of  umbral  oscillations,  then, 
we  can  look  for  free  eigenmodes  of  oscillation  in  an  isothermal 
layer  (layer  2)  bounded  above  and  below  by  perfectly  reflecting 
rigid  boundaries.  We  seek  solutions  to  the  wave  equations  (2), 
(3)  in  the  isothermal  region  0  <  z  <  h,  subject  to  the 
boundary  conditions 

u  =  w  =  0  at  z  =  0 ,  h  . 
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We  expect  that  the  calculated  eigenmodes  will  be  fairly 
independent  of  h  for  sufficiently  large  h  (h  >  10  Hq  for 
the  lowest  mode)  because  most  of  the  downward  reflection  will 
have  occurred  below  z  =  h. 

Our  method  of  solution  is  to  choose  a  value  of  h  and 
calculate  the  eigenmodes  and  eigenf requencies  using  our 

numerical  procedure.  The  calculations  are  then  repeated  for 
a  larger  value  of  h  to  assure  that  we  have  achieved  the 
desired  accuracy  (ft_  to  three  significant  figures) .  The 
resulting  eigenf  requencies  for  the  first  three  inodes  are 

listed  in  Table  1  for  comparison  with  the  resonant  frequencies 
S2r  determined  from  the  forced  solutions.  As  exoected,  the 
agreement  between  S2  and  S2  is  good.  The  computed  eigen- 
modes  resemble  very  closely  the  resonant  modes  shown  in 
Figure  3, 

An  even  simpler  analytical  approach  to  the  free  eigenmodes 
is  possible.  The  nature  of  the  eigenmodes  is  most  closely 
associated  with  the  horizontal  motions,  with  the  vertical 
motions  playing  a  passive  role.  This  suggests  that  we  can 
neglect  the  terms  involving  the  vertical  velocity  w  in  the 
equation  for  the  horizontal  velocity  u,  Equation  (2).  This 
gives  a  closed  equation  for  u,  and  the  vertical  motion  w 
is  obtained  by  solving  Equation  (3)  with  the  known  u(z). 

This  is  the  "quasi-Alfven"  approximation  of  Uchida  and 
Sakurai  (1975),  which  is  also  used  by  Hollweg  (1979)  .  These 
authors  justify  the  approximation  solely  on  the  basis  that 


2  2  . 

VA>>  c  9  w^lc^  1S  not  true  in  the  lower  part  of  the  atmosphere 

where  much  of  the  energy  of  the  eigenmodes  resides;  in  fact, 

2  2 

<  c  at  z  =  0  in  both  Uchida  and  Sakurai's  and  our  umbral 
model.  The  real  justification  for  the  approximation  low  in 
the  atmosphere  is  that  the  actual  computed  values  of  the  terms 
involving  w  in  Equation  (2)  are  much  smaller  than  terms 
involving  u  for  the  computed  eigenmodes.  High  in  the 

atmosphere  the  approximation  is  valid  on  the  basis  that 

2  ^  2 
VA>>C  • 


Under  the  quasi-Alfven  approximation,  Equation  (2)  reduces 


to 
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Ho 


(18) 


2  2 

in  an  isothermal  layer  where  v^  = v^Q  exp(z/HQ).  The  solution 
to  this  equation  satisfying  the  conditions  u(0)  =0  and 
u  +  0  as  z  *>  00  is  given  in  terms  of  Bessel  functions  in  the 
form 

"2c 
7Ao 


Gn(z)  =  J2K[v~(nn_I<2),5exp(-z/2Ho) 


(19) 


where 


a  = 

n 


Ao 


.2 


2  D2K,n  +  K 


(20) 


and  where  j_  is  the  nth  zero  of  J_  (z) .  We  consider 

z  ,  n  &  k 

the  eigenf requencies  ft  as  approximate  values  of  the 
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eigenfrequencies  computed  using  the  full  wave  equations 

and  denote  them  by  the  symbol  Values  of  for  the 

lowest  three  eigenmodes  are  given  in  Table  1 .  The  values  of 
nE  are  indeed  close  to  the  corresponding  values  of  and, 

coincidentally,  even  closer  to  the  corresponding  values  of 

Thus,  the  simple  quasi-Alfven  solution  (19) ,  (20) ,  which 
is  similar  to  that  of  Uchida  and  Sakurai  (the  upper  boundary 
condition  is  slightly  different) ,  gives  a  good  estimate  of 
the  oscillation  period  even  though  certain  details  of  the 
wave  motion  are  lost . 

6.  Conclusion  and  Discussion 

We  identify  umbral  oscillations  with  the  lowest  resonant 
mode  of  fast  magneto-atmospheric  wave  in  our  model  umbral 
atmosphere.  The  trapping  of  this  resonant  mode  is  due 
primarily  to  the  increasing  Alfven  speed  upward  into  the 
chromosphere  and  the  increase  in  sound  speed  downward  into 
the  convection  zone;  thus,  both  the  magnetic  and  acoustic 
nature  of  the  wave  are  quite  important,  with  buoyancy  playing 
a  smaller  role.  This  trapping  mechanism  is  similar  to  that 
proposed  for  running  penumbral  waves  by  Nye  and  Thomas  (1974, 
1976)  ,  for  the  case  of  a  horizontal  magnetic  field.  In  the 
present  case,  the  downward  reflection  is  not  complete;  a 
small  fraction  of  the  total  energy  in  the  mode  escapes  to 
large  heights  by  converting  into  the  form  of  an  acoustic 
wave  along  the  vertical  magnetic  field  lines.  This  small 


energy  leak  may  be  important  in  heating  the  corona  above  an 
active  region.  Downward  reflection  from  the  corona  is  unim¬ 
portant  in  determining  the  character  of  umbral  oscillations, 
since  most  of  their  energy  is  reflected  well  below  the  transi¬ 
tion  region. 

The  period  of  umbral  oscillation  in  our  model,  for  our 
assumed  values  of  the  parameters  (Section  2) ,  is  153  s.  If 
we  change  only  the  value  of  the  magnetic  field  strength,  we 
would  of  course  change  this  period.  However,  we  accept  the 
convincing  argument  of  Uchida  and  Sakurai  (1975)  that  over¬ 
stable  convecfion  occurs  only  at  a  level  where  the  Alfven 
speed  is  not  too  different  from  the  sound  speed.  Above  or 
below  this  level  the  magnetic  field  lines  are  too  stiff  or 
too  weak,  respectively,  for  overstable  oscillation  to  occur. 

This  means  that,  regardless  of  the  field  strength  of  the  umbra, 
the  value  of  the  Alfven  speed  v^Q  at  the  level  of  forcing  z  = 

0  is  roughly  the  same.  According  to  Uchida  and  Sakurai,  this 
accounts  for  the  fairly  narrow  range  of  observed  periods  of 
umbral  oscillation  among  spots  of  widely  differing  field  strength. 
In  terms  of  our  model,  it  means  that  we  represent  different 
sunspots  by  changing  the  field  strength  BQ  and  the  level  of 
forcing  z  =  0  simultaneously,  with  vAq,  and  hence  the  period  of 
oscillation,  nearly  unchanged.  The  observed  spread  of  oscillation 
periods  is  more  likely  due  to  variations  in  the  size  of  sunspots, 
that  is,  variations  in  the  horizontal  wavenumber  k. 
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Finally,  a  comment  should  be  made  concerning  the  possible 
cooling  of  sunspots  by  magnetohydrodynamic  waves.  One  of  us 
(Thomas  1978)  has  shown  that  pure  Alfven  waves  (k  =  0)  are 
reflected  strongly  downward  in  the  umbral  photosphere,  but 
are  only  weakly  reflected  upwards  in  the  umbral  convection  zone. 
This  left  open  the  possibility  of  cooling  by  downward  - 
propagating  waves.  However,  it  was  speculated  that  when  more 
realistic  magneto-atmospheric  waves,  including  the  effects  of 
compression,  are  considered,  the  upward  reflection  would  be 
much  stronger.  The  results  for  layer  1  in  the  present  paper 
confirm  this.  The  upward  reflection  of  waves  in  layer  1  is  very 
strong  for  the  horizontal  wavenumber  k  considered  here  and  also 
for  higher  values  of  k.  This,  together  with  the  results  of 
Thomas  (1978) ,  argues  strongly  against  significant  cooling  of 
sunspots  by  waves. 
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TABLE  1 


The  nondimensional  resonant  frequency  for  forced  oscillations 

(Section  4.1),  the  exact  free  eigenf requency  (Section  5),  and 

the  approximate  free  eigenfrequency  (Section  5) ,  for  the 

Jb 

lowest  three  modes  (n=l,2,3)  in  each  case.  The  dimensional 
periods  t  ,  corresponding  to  the  values  of  n_,  are  based  on 

K 

the  values  T  =  4500  K,  B  =  1000  G. 
o  o 


n 

xR  (s) 

1 

0.712 

0.743 

0.713 

153 

2 

1.51 

1.52 

1.50 

72 

3 

2.33 

2.46 

2.29 

47 

72 
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Fig.  1.  The  three-layer  model  of  the  umbral  atmosphere. 


Fig.  2 


Fig.  3 


Fig.  4 


Fig.  5 


Fig.  6 


Fig.  7 


.  Peaks  energy  density  E  as  a  function  of  the  forcing 

ITlaX 

frequency  ft  for  horizontal  forcing  in  the  two-layer 
model  atmosphere.  Each  plotted  point  represents  the 
result  of  a  numerical  solution.  The  numerical  values 
of  Emax  are  based  on  the  forcing  amplitude  fi(0)=l  km  s  . 

.  The  normalized  horizontal  velocity  as  a  function  of 
height  in  layer  2  for  the  lowest  three  resonant  fre¬ 
quencies  in  the  two-layer  model  atmosphere.  Each  higher 
mode  has  an  additional  node  and  all  the  profiles  have 
the  same  asymptotic  behavior  for  large  z. 

.  The  logarithmic  derivative  of  the  vertical  velocity 
as  a  function  of  height  in  layer  2  for  the  lowest  re¬ 
sonant  frequency  in  the  two- layer  model  atmosphere. 

.  The  total  energy  density  E  as  a  function  of  height  in 
layer  2  for  either  horizontal  or  vertical  forcing  at 
the  lowest  resonant  frequency  in  the  two-layer  model 
atmosphere. 

.  Horizontal  and  vertical  velocity  profiles  and  the  total 
energy  density  distribution  in  layer  1  for  horizontal 
forcing  (a,c)  and  vertical  forcing  (b,d)  at  z=0.  The 
horizontal  velocity  Q  is  the  solid  curve  and  the  verti¬ 
cal  velocity  Q  is  the  dashed  curve. 

.  The  total  energy  density  E  as  a  function  of  height  in 
layer  2  for  horizontal  (solid  curve)  and  vertical 
(dashed  curve)  forcing  at  ft=0.743,  about  4%  above  the 
lowest  resonant  frequency,  in  the  two-layer  model 


73 


atmosphere. 

Fig.  8.  The  total  energy  density  E  as  a  function  of  height  in 

layers  2  and  3  for  either  horizontal  or  vertical 

forcing  in  the  three-layer  model  atmosphere.  The  energy 

density  drops  three  orders  of  magnitude  across  the 

transition  region  z=z  =20  H  and  is  not  visible  above 

t  o 

z=zfc  at  this  scale. 


